\ 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