# rkf45_sys.f

\ 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,

[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
;

\ 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]