# rkf45a.f

\ 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])"
;