\ Numerical solution of first-order ordinary differential equation \ by 4th order Milne predictor-corrector algorithm following \ Abramowitz & Stegun p. 896, 25.5.13 \ Solves dx/dt = f(x,t) \ P: x(n+1) = x(n-3) + (4h/3) * [2f(n) - f(n-1) +2f(n-2)] \ C: x(n+1) = x(n-1) + (h/3) * [f(n-1) + 4f(n) + f(n+1)] \ This is an ANS Forth program requiring FLOATING and FLOATING EXT \ wordsets. This program assumes a separate floating point stack. \ For simplicity and clarity, the program is intended to be compiled \ using a FORmula TRANslator \ Example \ : fnb ( f: x t -- t^2*exp[-x]) \ f^2 FSWAP FNEGATE FEXP F* ; \ use( fnb 0e0 0e0 5e0 0.1e0 )milne MARKER -milne INCLUDE ftran111.f : fvariables 0 DO FVARIABLE LOOP ; 5 fvariables t h h3 h43 tf 8 fvariables x x0 xm1 xm2 xm3 y y0 ym1 3 fvariables f0 fm1 fm2 4 fvariables xk1 xk2 xk3 xk4 v: fdummy : update \ shuffle variables f" xm3 = xm2" f" xm2 = xm1" f" fm2 = fm1" f" xm1 = x0" f" ym1 = y0" f" fm1 = f0" f" x0 = x" f" y0 = y" ; : runge_step f" f0 = fdummy(x,t)" f" xk1 = h*f0/2 " \ compute k1 f" t = t +h/2 " \ t -> t + h/2 f" xk2 = h*fdummy(x+xk1,t)" f" xk3 = h*fdummy(x+xk2/2,t)" f" t = t+h/2 " \ t -> t + h/2 f" xk4 = h*fdummy(x+xk3,t)/2" f" x = x + (xk1 + xk2 + xk3 + xk4)/3 " ; : x' \ integration step f" f0 = fdummy(x,t)" f" x = xm3 + h43 * (2*(f0 + fm2) - fm1)" f" t = t + h" f" y = ym1 + h3 * (fm1 + 4* f0 + fdummy(x,t))" ; : display CR t F@ FS. 5 SPACES x F@ FS. 5 SPACES f" ln(1+t^3/3)" FS. ; : done? t F@ tf F@ F> ; : initialize ( xt --) ( f: x0 t0 tf h -- ) defines fdummy \ vector fn_name h F! tf F! t F! x F! \ initialize variables f" h3 = h/3" f" h43 = 4*h3" 0.e0 FDUP xm2 F! \ put 0 in these FDUP xm1 F! FDUP x0 F! FDUP y0 F! fm1 F! 4 0 DO runge_step f" y = x" update LOOP ; : )milne ( xt --) ( f: x0 t0 tf h -- ) initialize BEGIN display x' \ indefinite loop done? UNTIL ;