\ Numerical solution of a system of first-order ordinary
\ differential equations 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 are n-dimensional
\ vectors.
\ For simplicity and clarity, the program is intended to be
\ compiled using a FORmula TRANslator
FALSE [IF]
Example:
: fnb ( f: t -- ) ( [f] [x] --)
LOCALS| x{ f{ |
FRAME| aa |
x{ _len 0 DO
f" f{ I } = aa^2 * exp(-x{ I })"
LOOP
|FRAME
;
x 3 tab->vec 0e0, 0e0, 0e0, \ initial values
use( fnb 0e0 5e0 0.1e0 eps: 1e-5 )rkf45_sys
[THEN]
FALSE [IF]
Example:
: fnb ( f: t -- ) ( [f] [x] --)
LOCALS| x{ f{ |
FDROP \ no t dependence
f" f{ 0 } = x{ 1 }"
f" f{ 1 } = -x{ 0 }"
f" f{ 2 } = 0"
;
x 3 tab->vec 0e0, 1e0, 0e0, \ initial values
use( fnb 0e0 5e0 0.1e0 eps: 1e-5 )rkf45_sys
[THEN]
MARKER -runge
include ftran202.f
\ uses a FORmula TRANslator for clarity
include arrays.f
\ uses "Scientific Forth" array notation
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 and FVALUEs
\ can be passed to expressions by enclosing them in [ ].
\ ---------------------------- 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
: fvariables: 0 DO FVARIABLE LOOP ;
6 fvariables: t tf epsilon delta MaxErr h
3 VALUE Neq
Neq long 1 FLOATS 1array x
Neq long 1 FLOATS 1array y
Neq long 1 FLOATS 1array k0
Neq long 1 FLOATS 1array k1
Neq long 1 FLOATS 1array k2
Neq long 1 FLOATS 1array k3
Neq long 1 FLOATS 1array k4
Neq long 1 FLOATS 1array k5
Neq long 1 FLOATS 1array f \ holds function values
FALSE [IF]
Note: I violate my vector naming convention, dropping the
final "{", because I will retain it in LOCAL names
within the functions.
[THEN]
v: fdummy
\ -------------------------------------------- debugging code
0 VALUE debugging
\ TRUE TO debugging
: debug\ debugging 0= IF POSTPONE \ THEN ; IMMEDIATE
\ ---------------------------------------- end debugging code
FALSE [IF]
As a courtesy I have provided error checking code via `debug\'
above. If the VALUE `debugging' has been set to TRUE, the line
following `debug\' will be compiled. Otherwise it will not.
It is the programmer's responsibility to enforce type- and length
compatibility in the vector functions below (scalar multiplication,
vector addition).
[THEN]
: eps: ; \ syntactic sugar
: compatible ( [v1] [v2] --)
_len SWAP _len <> ABORT" incompatible vectors"
;
: Smul ( [lambda] [v1] [v2] ) \ scalar multiplication
\ v2 = lambda * v1
LOCALS| v2{ v1{ lambda |
debug\ v1{ v2{ compatible
v1{ _len 0 DO f" v2{ I } = lambda * v1{ I }" LOOP
;
: Vadd ( [v1] [v2] [v3] -- ) \ vector addition
\ v3 = v1 + v2
LOCALS| v3{ v2{ v1{ |
debug\ v1{ v2{ compatible v2{ v3{ compatible
v3{ _len 0 DO f" v3{ I } = v1{ I } + v2{ I }" LOOP
;
: Lcomb ( [v1] [v2] [v3] -- ) ( f: lambda mu --)
\ linear combination: v3 = lambda * v1 + mu * v2
LOCALS| v3{ v2{ v1{ |
FRAME| aa bb |
debug\ v1{ v2{ compatible v2{ v3{ compatible
v3{ _len 0 DO f" v3{ I } = bb*v1{ I } + aa*v2{ I }" LOOP
|FRAME
;
: display ( [x] )
f" t" F. 5 SPACES
LOCALS| x{ |
x{ _len 0 DO f" x{ I }" FS. 3 SPACES LOOP
;
: do_k0 ( [x] [k0] ) \ k0 = h*f(x,t)
f" fdummy([f],[x],t)"
f" Smul([h],[f],[k0])"
debug\ CR ." k0: " f" display([k0])"
;
: do_k1 \ k1 = h*f(x+k0/4,t+h/4)
1e0 a10 f" Lcomb( [x], [k0], [y])" \ y = 1*x + a10 * k0
f" fdummy( [f], [y], t+[c1]*h )" \ f = dummy(y,t+c1*h)
f" Smul([h],[f],[k1])" \ k1 = h*f
debug\ CR ." k1: " f" display([k1])"
;
: do_k2 \ k2 = h*f( x + a20*k0 + a21*k1, t + 3*h/8 )
a21 a20 f" Lcomb( [k1], [k0], [y] )" \ y = a20*k0 + a21*k1
f" Vadd( [x], [y], [y] )" \ y = x + y
f" fdummy([f],[y],t+[c2]*h)" \ f = dummy(y,t+c2*h)
f" Smul([h],[f],[k2])" \ k2 = h* f
debug\ CR ." k2: " f" display([k2])"
;
: do_k3 \ k3 = h*f( x + a30*k0 + a31*k1 + a32*k2, t + 12*h/13 )
x y k0 k1 k2 LOCALS| k2{ k1{ k0{ y{ x{ |
x{ _len 0 DO
f" y{ I } = x{ I } + [a30]*k0{ I } + [a31]*k1{ I } + [a32]* k2{ I }"
LOOP
f" fdummy([f],[y],t+[c3]*h)"
f" Smul([h],[f],[k3])"
debug\ CR ." k3: " f" display([k3])"
;
: do_k4 \ k4 = h*f( x + a40*k0 + a41*k1 + a42*k2 + a43*k3, t + h)
x y k0 k1 k2 k3 LOCALS| k3{ k2{ k1{ k0{ y{ x{ |
x _len 0 DO
f" y{ I } = x{ I } + [a40]*k0{ I } + [a41]*k1{ I }
+ [a42]*k2{ I } + [a43]*k3{ I }"
LOOP
f" fdummy([f],[y],t+h)"
f" Smul([h],[f],[k4])"
debug\ CR ." k4: " f" display([k4])"
;
: do_k5 \ k5 = h*f( x + a50*k0 + a51*k1 + a52*k2 +
\ a53*k3 + a54*k4, t + h/2)
x y k0 k1 k2 k3 k4 LOCALS| k4{ k3{ k2{ k1{ k0{ y{ x{ |
x{ _len 0 DO
f" y{ I } = x{ I } + [a50]*k0{ I } + [a51]*k1{ I }
+ [a52]*k2{ I } + [a53]*k3{ I } + [a54]*k4{ I }"
LOOP
f" fdummy([f],[y],t+0.5*h)"
f" Smul([h],[f],[k5])"
debug\ CR ." k5: " f" display([k5])"
;
: DeltaX ( [k0] [k2] [k3] [k4] [k5] --) ( f: -- DeltaX)
LOCALS| k5{ k4{ k3{ k2{ k0{ |
f" delta=0"
k0{ _len 0 DO
f" delta = delta
+ ( [b0]*k0{ I } + [b2]*k2{ I }
+ [b3]*k3{ I } + [b4]*k4{ I }
+ [b5]*k5{ I } )^2 "
LOOP
f" sqrt( delta )"
;
: MoreThan POSTPONE F> ; IMMEDIATE
: step_too_big?
f" delta = DeltaX( [k0], [k2], [k3], [k4], [k5])"
f" MaxErr = abs(h)*epsilon"
f" MoreThan(delta,MaxErr)"
;
: 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 ( [x] [k0] [k1] [k2] [k3] [k4] [k5] --)
LOCALS| k5{ k4{ k3{ k2{ k1{ k0{ x{ |
x{ _len 0 DO
f" x{ I } = x{ I } + [d0]*k0{ I }+[d2]*k2{ I }
+[d3]*k3{ I }+[d4]*k4{ I }+[d5]*k5{ I }"
LOOP
f" t = t+h"
;
: step \ integration step
do_k0 do_k1 do_k2 do_k3 do_k4 do_k5
step_too_big?
IF \ error too large, reduce step size and repeat
shrink \ reduce h
\ however, don't make h' < 0.25*h
RECURSE \ repeat step with smaller |h|
ELSE \ error acceptable
f" update([x], [k0], [k1], [k2], [k3], [k4], [k5])"
\ error too small? perhaps increase step size
expand
THEN
;
: initialize ( xt --) ( f: t0 tf h epsilon -- )
f" epsilon=" f" h=" f" tf=" f" t="
defines fdummy \ vector fn_name
;
: display
CR t F@ FS.
x _len 0 DO
3 SPACES x I } F@ FS.
LOOP
\ 3 SPACES f" ln(1+t^3/3)" FS.
3 SPACES f" sin(t)" FS.
;
: not_done? t F@ h F@ F+ tf F@ F< ;
: )rkf45_sys ( xt --) ( f: t0 tf h -- )
initialize
BEGIN not_done? \ solve
WHILE step display
REPEAT
f" h = tf - t"
step
display
;
\ initialization
: get_f# ( f: -- x)
[CHAR] , PARSE >FLOAT 0= ABORT" Wrong fp# format!" ;
: TAB->VEC ( V{ n -- ) \ input vector from a table
LOCALS| N v{ |
N 0 DO get_f# v{ I } F! LOOP
;
FALSE [IF]
\ Example:
\ f{ is the vector of function values
3 fvariables: sigma r beta
: fnb ( [f] [x] -- ) ( f: t -- )
LOCALS| x{ f{ |
FDROP
f" f{ 0 } = sigma * ( x{ 1 } - x{ 0 })"
f" f{ 1 } = (r - x{ 2 }) * x{ 0 } - x{ 1 }"
f" f{ 2 } = x{ 0 } * x{ 1 } - beta * x{ 2 }"
;
x 3 tab->vec 2e0, 5e0, 5e0, \ initial values for Lorenz model
f$" sigma = 10" f$" beta = 8/3" f$" r = 28"
\ use( fnb 0e0 1e0 0.1e0 eps: 1.e-5 )rkf45_sys
[THEN]