Rung-Kutta-Fehlberg (RKF) method is a numerical method of solving ordinary differential equations derived from the Runge-Kutta method. RKF uses a fourth order and fifth order accurate method to provide 2 estimates of the next step in the solution; the difference in these results is related to the accuracy of the solution and can be used to decrease the integration step length if the error exceeds a specified accuracy. The clever feature of RKF is that it reduces the number of evaluations of the force functions by using an integrator that uses some of the same values in both the 4th order and 5th order methods.
The attached code compares RKF solution with the TCLlib::math calculus RK4 4th order method. The RKF implementation is designed to emulate the RK4 code with the addition of the error vector to specify the desired accuracy.
gwm 17-05-05 correction to recommended step length improves speed, and accuracy. GWM 25-05-05 add buttons to exercise the solver, and remove code common to Another Graphing Widget. The plot code from that page should be includes as a separate file which I call plotcanvas.tcl, but you can either - include the code directly to replace the source plotcanvas.tcl or save as a tcl file.
plotcanvas.tcl is the contents of the page Another Graphing Widget which can be found in mini.net/tcl and is used for display of results. The code is not duplicated here.
source plotcanvas.tcl
# rk-f properties, see http://math.fullerton.edu/mathews/n2003/RungeKuttaFehlbergMod.html
# basicaly replaces 4th order with 4th and fifth order solution and
# calculates estimate of time step required for next step of integration.
# By careful selection of the RK coefficients, fewer evaluations of derivatives etc are required.
proc rungeKuttaFehlbergStep { t tstep xvec func {errvec ""}} {
# find steps required to reach tolerance errvec (for each element of the xvec state function)
# divide step tstep into N equally spaced steps (tstep/N) if estimated step length < tstep.
#
# Six steps:
# - k1=tstep*f(t,y)
# - k2=tstep*f(t+.25.tstep,y+.25*k1)
# - k3=tstep*f(t+3tstep/8, y+(3.k1+9.k2)/32) => (tt,y+3/8*(k1/4+3k2/4))
# - k4=tstep*f(t+12tstep/13, y+(1932.k1-7200.k2+7296*k3)/2197)
# - k5=tstep*f(t+tstep, y+(439.k1/216-8.k2+3680.k3/513-845*k4/4104))
# - k6=tstep*f(t+tstep/2, y+(-8.k1/27+2.k2-3544.k3/2565+1859.k4/4104-11.k5/40))
# 4th order sol is y4th=y+25/216*k1+1408/2565.k3+2197/4104*k4-k5/5
# 5th order sol is y5th=y+16/135*k1+6656/12825.k3+28561/56430*k4-9k5/50+2/55k6
# and opt step size is s.h where err is the desired accuracy (dimension cancels with y5,y4)
# s = {err.h/(2*(y5th-y4th)))}^0.25
# or [s^4*2*(y5th-y4th)]/h=err??
#
set funcq [uplevel 1 namespace which -command $func]
if {$errvec != ""} {
upvar 1 $errvec errv ;# set "local" address of variable with constant error desired.
} else {set errv $xvec}
# - k1=f(t,y)
set xk1 [$funcq $t $xvec] ;# first step is vector of derivatives f(t,y)
set tstep2 [expr {$tstep/4.0}]
set xvec2 {}
foreach x1 $xvec xv1 $xk1 {
lappend xvec2 [expr {$x1+$tstep2*$xv1}]
}
# - k2=tstep*f(t+.25.tstep,y+.25*k1) tstep2=dt/4
set xk2 [$funcq [expr {$t+$tstep2}] $xvec2]
# - k3=tstep*f(t+3tstep/8, y+(3.k1+9.k2)/32) => xk3=tstep*f(t+tstep3,y+3/8*(k1/4+3k2/4))
# =tstep3*f(t+tstep3,y+(k1/4+3k2/4))
set tstep3 [expr {3.0*$tstep/8.0}]
set xvec3 {}
foreach x1 $xvec xv1 $xk1 xv2 $xk2 {
lappend xvec3 [expr {$x1+$tstep3*(.25*$xv1+.75*$xv2)}]
}
set xk3 [$funcq [expr {$t+$tstep3}] $xvec3]
# - k4=tstep*f(t+12tstep/13, y+(1932.k1-7200.k2+7296*k3)/2197) =>
# xk4=tstep4*f(t+tstep4, y+(161.k1-600.k2+608*k3)/169)
set tstep4 [expr {12.0*$tstep/13.0}]
set xvec4 {}
foreach x1 $xvec xv1 $xk1 xv2 $xk2 xv3 $xk3 {
lappend xvec4 [expr {$x1+$tstep4*(161.0*$xv1-600.0*$xv2+608.0*$xv3)/169.0}]
}
set xk4 [$funcq [expr {$t+$tstep4}] $xvec4]
# k5=tstep*f(t+tstep, y+(439.k1/216-8.k2+3680.k3/513-845.0*k4/4104)) =>
# xk5=tstep5*f(t+tstep5, y+(439.k1/216-8.k2+3680.k3/513-845*k4/4104))
set tstep5 $tstep
set xvec5 {}
foreach x1 $xvec xv1 $xk1 xv2 $xk2 xv3 $xk3 xv4 $xk4 {
lappend xvec5 [expr {$x1+$tstep5*(439./216*$xv1-8.*$xv2+3680.0/513*$xv3-845.0/4104*$xv4)}]
}
set xk5 [$funcq [expr {$t+$tstep5}] $xvec5]
# - k6=tstep*f(t+tstep/2, y+(-8.k1/27+2.k2-3544.k3/2565+1859.k4/4104-11.k5/40)) =>
# xk6=tstep6*f(t+tstep/2, y+(-16.k1/27+2.k2-7088.k3/2565+1859.k4/2052-11.k5/20))
set tstep6 [expr {$tstep/2.0}]
set xvec6 {}
foreach x1 $xvec xv1 $xk1 xv2 $xk2 xv3 $xk3 xv4 $xk4 xv5 $xk5 {
lappend xvec6 [expr {$x1+$tstep*(-8./27*$xv1+2.*$xv2-3544./2565*$xv3+1859./4104.*$xv4-11./40*$xv5)}]
}
set xk6 [$funcq [expr {$t+$tstep6}] $xvec6]
# 4th order sol is y4th=y+25/216*k1+1408/2565.k3+2197/4104*k4-k5/5
# 5th order sol is y5th=y+16/135*k1+6656/12825.k3+28561/56430*k4-9k5/50+2/55k6
# use 5th order result
set result {}
set res4 {}
set stepest $tstep
foreach x0 $xvec k1 $xk1 k3 $xk3 k4 $xk4 k5 $xk5 k6 $xk6 err $errv { ;# values {k2 $xk2} are not used
set dx [expr {16.0/135.0*$k1+6656.0/12825.0*$k3+28561.0/56430.0*$k4-0.18*$k5+2.0/55.0*$k6}]
lappend result [expr {$x0+$dx*$tstep}]
if {$errvec != ""} {
set dx4 [expr {25.0/216.0*$k1+1408.0/2565*$k3+2197/4104.0*$k4-0.2*$k5}]
lappend res4 [expr {$x0+$dx4*$tstep}]
if {$dx!=$dx4} {
#step size recommended is s = {err.h/(2*(y5th-y4th)))}^0.25
set spp [expr {$tstep*pow($err*$tstep/(2.0*$tstep*(abs( $dx-$dx4 ))) , 0.25)}]
if {$spp <$stepest} {
set stepest $spp
}
}
}
}
if {$stepest<$tstep} { ;# adjust time step by redoing N sub-steps
set nsubsteps 2;#[expr {int(ceil( $tstep/$stepest))}] ;# or 2
set dtsub [expr {double($tstep)/$nsubsteps}]
set result $xvec ;# restart half time step
set i 0
set tsu $t
while {$i<$nsubsteps} {
set result [rungeKuttaFehlbergStep $tsu $dtsub $result $func errv ]
# if the estimate is not good enough then we will recurse again in the sub-step
set tsu [expr {$tsu+$dtsub} ]
incr i
}
}
return $result
}if {1==0} { These procs exercise the RKF solver
red line is the RKF tme step adjusted solution (there may be more than 1 step between each graph point) green line is the RK4 solution using the base timesteps of the RKF solution
and compares to a short step RK4 solver (curve in black).
The short step RK4 is assumed to be accurate. The RKF points coincide with the short step RK4 curve. The long step RK4 points do not coincide. The RKF algorithm has automatically adjusted the time steps to (1) fit exactly N steps between the basic steps and (2) keep the estimated error below the requested accuracy.
This technique is useful for dynamics calculations where long time steps are needed for the motion between collisions, and short time steps during a collision. Also for calculating the trajectory of a highly eccentric object (such as a comet) around the sun where shorter steps are useful when the object is near the sun as the acceleration of the object is high. (EG Halley's comet spends 75 years orbiting the sun, of which less than 3 months is typically visible from earth.)
}
proc damposc {} {
# test with a damped oscillator; long time steps result in error using RK-4th order,
# RKF automatically halves steps to retain a maximum error.
package require math::calculus ;# used for comparison with supplied RK4 algorithm
proc dampened_oscillator {t xvec} { ;# as used by RungeKuttaStep also used by RKFstep.
global damp
set x [lindex $xvec 0] ;# pos
set x1 [lindex $xvec 1] ;# velocity
# dpos/dx=v
# dv/dx= k.(x-origin)-damp*vel
return [list $x1 [expr {-$damp*$x1-$x}]]
}
global damp
set tmax 40.0
set damp 0.0
catch {destroy .fib}
set dr2 [grapharea .fib -width 600 -height 400]
$dr2 configure -pencolour gray
$dr2 configure -xmin -1 -xmax $tmax -ymin -1.2 -ymax 1.2
pack $dr2
$dr2 configure -pencolour gray50
$dr2 showgrid 5 1
$dr2 configure -pencolour orange
while {$damp<1} {
puts "Damp is $damp"
set xvec { 1.0 0.0 } ;# position, velocity
set xvec4 { 1.0 0.0 } ;# position, velocity rk4 solver
set t 0.0
set errvec {0.001 0.0002}
set tstep 2.5 ;# if tstep <1.5 then RK4 with constant time step is acceptable
# try tstep 3 to see rk4 diverging drastically!
set pos ""
lappend pos $t [lindex $xvec 0]
set pos4rk ""; lappend pos4rk $t [lindex $xvec4 0]
while { $t < $tmax } {
set res4 [::math::calculus::rungeKuttaStep $t $tstep $xvec4 dampened_oscillator]
set result [rungeKuttaFehlbergStep $t $tstep $xvec dampened_oscillator errvec]
set t2 [expr {$t+$tstep}]
lappend pos $t2 [lindex $result 0]
lappend pos4rk $t2 [lindex $res4 0]
set t $t2
set xvec $result
set xvec4 $res4
}
$dr2 configure -pencolour red
$dr2 create line $pos
$dr2 create text 10.0 1.1 -text {{Runge-Kutta-Fehlberg error correcting}}
$dr2 configure -pencolour green
$dr2 create line $pos4rk
$dr2 create text 10.0 0.9 -text { "Runge-Kutta long step"}
# redo calculation using RK 4th order, but very short time step
set xvec4 { 1.0 0.0 } ;# position, velocity rk4 solver
set t 0
set tstep [expr {$tstep/10.0}]
set pos "";lappend pos $t [lindex $xvec4 0]
while { $t < $tmax } {
set res4 [::math::calculus::rungeKuttaStep $t $tstep $xvec4 dampened_oscillator]
set t2 [expr {$t+$tstep}]
lappend pos $t2 [lindex $res4 0]
set t $t2
set xvec4 $res4
}
$dr2 configure -pencolour black
$dr2 create line $pos
$dr2 create text 10.0 1.0 -text {"Runge-Kutta short step"}
set damp [expr {$damp+0.1}]
update
}
}
proc solidwall {} {
puts "Solidwall test"
# hitwal - 2 solid walls with ball bouncing between.
proc wallhit {t xvec} { ;# as used by RungeKuttaStep also used by RKFstep.
set x [lindex $xvec 0] ;# pos
set x1 [lindex $xvec 1] ;# velocity
# dpos/dx=v
# dv/dx= if x>wall k.(x-wall) or if x<-wall f=-k.wall-x
set wc 1.e10
set f 0
if {$x>1.0} { set f [expr {-$wc*($x-1.0)}]
} elseif {$x<-1.0} { set f [expr {-$wc*($x+1.0)}]
}
#puts "$t $xvec -> forces $x1 $f"
return [list $x1 $f]
}
set xvec { 0.0 0.5 } ;# position, velocity
set errvec {0.0002 2.e-5}
set t 0.0
set tmax 22.2
catch {destroy .fib}
set dr2 [grapharea .fib -width 600 -height 400]
$dr2 configure -pencolour gray
$dr2 configure -xmin -1 -xmax $tmax -ymin -1.2 -ymax 1.2
pack $dr2
$dr2 configure -pencolour gray50
$dr2 showgrid 5 1
$dr2 configure -pencolour orange
set tstep 0.01 ;# take 10 steps to hit wall
set pos ""
while {$t<$tmax} {
set result [rungeKuttaFehlbergStep $t $tstep $xvec wallhit errvec]
set xvec $result
puts "Time $t [lindex $result 0] [lindex $result 1]"
lappend pos $t [lindex $result 0]
set t [expr {$t+$tstep}]
update
}
$dr2 create line $pos
}
# create buttons for demos
set entrypts {}
puts "run solidwall to exercise the RKF collision with solid wall test"
puts "run damposc to see solutions to damped oscillator"
lappend entrypts {solidwall "exercise the RKF collision with solid wall test"}
lappend entrypts {damposc "solutions to damped oscillator"}
catch {destroy .testrkf}
set fex [frame .testrkf]
foreach ep $entrypts {
set choice [lindex $ep 0]
button $fex.$choice -text "[lindex $ep 1]" -command [lindex $ep 0]
pack $fex.$choice -side left
}
pack $fex -side topSee the reference [L1 ] for more details of RKF method.