\ Numerical solution of a first-order ordinary
\ differential equation by 4th+5th -order Runge-Kutta
\ -Fehlberg algorithm, following Paul L. DeVries, "A First
\ Course in Computational Physics" (Wiley, New York, 1994),
\ p. 220ff.
\ ---------------------------------------------------
\ (c) Copyright 2005 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 dx/dt = f(x,t) where x and f scalars.
\ For simplicity and clarity, the program is intended to be
\ compiled using a FORmula TRANslator
FALSE [IF]
Example
: fnb ( f: x t -- t^2*exp[-x])
FRAME| aa bb |
f" aa^2 * exp(-bb)"
|FRAME
;
f$" x=0" use( fnb 0e0 5e0 0.1e0 eps: 1e-5 )runge45
[THEN]
MARKER -runge
INCLUDE ftran202.f
\ uses a FORmula TRANslator for clarity
INCLUDE flocals.f
\ See the file FTRANDOC.TXT for the deviations from the usual
\ FORTRAN usages. In particular, a variable name enclosed in
\ square brackets [ ] means "pass-by-reference" whereas one
\ without [ ] means "pass-by-value". FCONSTANTs are "pass-by-
\ reference".
: eps: ; \ syntactic sugar
: fvariables: 0 DO FVARIABLE LOOP ;
6 fvariables: t tf epsilon delta MaxErr h
7 fvariables: x k0 k1 k2 k3 k4 k5
\ ---------------------------- define constants
f$" 1/4" FDUP FCONSTANT c1
FCONSTANT a10
f$" 3/8" FCONSTANT c2
f$" 3/32" FCONSTANT a20
f$" 9/32" FCONSTANT a21
f$" 12/13" FCONSTANT c3
f$" 1932/2197" FCONSTANT a30
f$" -7200/2197" FCONSTANT a31
f$" 7296/2197" FCONSTANT a32
f$" 439/216" FCONSTANT a40
-8e0 FCONSTANT a41
f$" 3680/513" FCONSTANT a42
f$" -845/4104" FCONSTANT a43
f$" 1/2" FCONSTANT c5
f$" -8/27" FCONSTANT a50
2e0 FCONSTANT a51
f$" -3544/2565" FCONSTANT a52
f$" 1859/4104" FCONSTANT a53
f$" -11/40" FCONSTANT a54
f$" 1/360" FCONSTANT b0
f$" -128/4275" FCONSTANT b2
f$" -2197/75240" FCONSTANT b3
f$" 1/50" FCONSTANT b4
f$" 2/55" FCONSTANT b5
f$" 16/135" FCONSTANT d0
f$" 6656/12825" FCONSTANT d2
f$" 28561/56430" FCONSTANT d3
f$" -9/50" FCONSTANT d4
f$" 2/55" FCONSTANT d5
\ -------------------- end constant definitions
v: fdummy
: MoreThan POSTPONE F> ; IMMEDIATE
: step_too_big?
f" delta = abs( [b0]*k0+[b2]*k2+[b3]*k3+[b4]*k4+[b5]*k5 )"
f" MaxErr = abs(h)*epsilon"
f" MoreThan(delta,MaxErr)"
;
: do_k0 f" k0 = h*fdummy(x,t)" ;
: do_k1 f" k1 = h*fdummy(x+[a10]*k0,t+[c1]*h)" ;
: do_k2 f" k2 = h*fdummy(x+[a20]*k0+[a21]*k1,t+[c2]*h)" ;
: do_k3 f" k3 = h*fdummy(x+[a30]*k0+[a31]*k1+[a32]*k2,t+[c3]*h)" ;
: do_k4 f" k4 = h*fdummy(x+[a40]*k0+[a41]*k1+[a42]*k2+[a43]*k3,t+h)" ;
: do_k5 f" k5 = h*fdummy(x+[a50]*k0+[a51]*k1+[a52]*k2+[a53]*k3+
[a54]*k4,t+[c5]*h)" ;
: shrink \ reduce step size
f" h = max( 0.9*h*sqrt(sqrt(Maxerr/delta)) , [c1]*h)"
\ however, don't make h' < 0.25*h
;
: expand \ increase step size
f" h = min( 0.9*h*sqrt(sqrt(Maxerr/delta)) , 4*h)"
\ however, don't make h' > 4*h
;
: update
f" x = x + [d0]*k0+[d2]*k2+[d3]*k3+[d4]*k4+[d5]*k5"
f" t = t+h"
;
: step \ integration step
do_k0 do_k1 do_k2 do_k3 do_k4 do_k5
step_too_big?
IF shrink \ reduce |h|
RECURSE \ repeat step with smaller |h|
ELSE \ error acceptable, update variables
update
\ error too small? perhaps increase step size
expand
THEN
;
: LessThan POSTPONE F< ; IMMEDIATE
: initialize ( xt --) ( f: t0 tf h epsilon -- )
8 set-precision
f" epsilon=" f" h=" f" tf=" f" t="
defines fdummy \ vector fn_name
;
: display ( [x] --)
LOCALS| x |
CR t F@ FS. 5 SPACES x F@ FS.
5 SPACES f" ln(1+t^3/3)" FS.
;
: not_done? t F@ h F@ F+ tf F@ F< ;
: )runge45 ( xt --) ( f: t0 tf h -- )
initialize
BEGIN not_done? \ solve
WHILE step f" display([x])"
REPEAT
f" h = tf - t"
step
f" display([x])"
;