\ Numerical solution of second-order linear differential equation
\ by (6th order) Numerov algorithm following
\ Hamming, p. 215.


\ ---------------------------------------------------
\     (c) Copyright 1999  Julian V. Noble.          \
\       Permission is granted by the author to      \
\       use this software for any application pro-  \
\       vided this copyright notice is preserved.   \
\ ---------------------------------------------------

\ This is an ANS Forth program requiring the
\   FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.

\ Environmental dependencies:
\       Assumes independent floating point stack


\ Solves      y'' = f(x)y + g(x)
\ via       Y_n+1 = 2*Y_n - Y_n-1 +
\                 + h^2 * [ f_n * Y_n + (g_n+1 - 2*g_n + g_n-1 )/12 ]
\           y_n = Y_n / (1 - h^2 * f_n /12 )

\    For simplicity and clarity, the program is intended to be compiled
\    using a FORmula TRANslator

\ Usage     use( f1 & f2 y y1 x0 xf h )numerov

\ Example
\    : fna   ( f: x -- -1.0 )   FDROP  -1e0  ;
\    : gna   ( f: x -- 0.0)     FDROP   0e0  ;
\    use( fna & gna 0e0  0.0999166111607e0 0e0 5e0 0.1e0 )numerov


MARKER  -numerov

INCLUDE ftran111.f

: &    '  ;     ( -- xt)

: fvariables     0 DO  FVARIABLE  LOOP  ;
13 fvariables h xf x y y0 y1 f g g0 g1 hsq hsq12 yy


v: fdummy
v: gdummy

: update    \ shuffle variables
    f" y0 = y"  f" g0 = g"
    f" y = y1"  f" g = g1"
;

: y'    \ integration step
    f" f  = fdummy(x)"
    f" x  = x+h"
    f" g1 = gdummy(x)"
    f" yy = y/(1 - hsq12 * f)"
    f" y1 = 2*y - y0 + hsq*f*yy + hsq12 * (g1 - 2*g + g0)"
;


: display       CR   f" x"  FS.
                5 SPACES   yy F@   FS.
                5 SPACES   f" sin(x)"  FS.
;


: done?  x F@  xf F@  F>  ;


: )numerov   ( xt --)  ( f: x0 t0 tf h -- )
        defines gdummy              \ vector fn_names
        defines fdummy
        h F!  xf F!  x F!           \ initialize variables
        y F!  y0 F!
        \ f" x = x+h"
        f" hsq = h^2"
        f" hsq12 = hsq / 12"
        BEGIN   y'  display         \ indefinite loop
                update
        done?   UNTIL  ;


\ Example: simple harmonic oscillator without driving term
: fna   ( f: x -- -1.0 )   FDROP  -1e0  ;
: gna   ( f: x -- 0.0)     FDROP   0e0  ;
use( fna & gna 0e0 0.0999166111607e0 0e0 5e0 0.1e0 )numerov