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








  HTMLized by Forth2HTML