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 implemetnation is designed to emulate the RK4 code with the addition of the error vector to specify the desired accuracy.
# rk-f properties, see http://math.fullerton.edu/mathews/n2003/RungeKuttaFehlbergMod.html
# basically 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
# 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} { ;# zero error gives an infinite best step size.
#step size recommended is s = {err.h/(2*(y5th-y4th)))}^0.25
set spp [expr {pow($err*$tstep/(2.0*(abs( $dx-$dx4 ))) , 0.25)}]
if {$spp<$stepest} {
set stepest $spp
}
}
}
}
if {$stepest<$tstep} { ;# adjust time step by redoing N sub-steps
set nsubsteps [expr {int(ceil( $tstep/$stepest))}]
set dtsub [expr double($tstep)/$nsubsteps]
# puts "at $t dt is $tstep rec step $stepest - use $nsubsteps steps $dtsub"
set result $xvec
set i 0
while {$i<$nsubsteps} {
set result [rungeKuttaFehlbergStep $t $dtsub $result $func errv ]
# if the estimate is not good enough then we will recurse again in the sub-step
set t [expr $t+$dtsub]
incr i
}
}
return $result
}
if {1==0} {This section is the graph drawing widget option found under Another Graphing Widget
}
proc grapharea {w args} { ;# a graph plotter canvas derived thing.
set dx 400
set dy 400
global $w.pencolour $w.xmin $w.xmax $w.ymin $w.ymax
set $w.pencolour black ;# as in Graphing Widget - intialise range of xy correctly
set $w.xmin -1; set $w.xmax 1
set $w.ymin -1; set $w.ymax 1
array set options {{-colour} red}
# split off the custom options:
set textArgs [list]
foreach {opt val} $args {
switch -- $opt {
{-pencolour} {set options($opt) $val}
default {lappend textArgs $opt $val}
}
}
eval canvas $w $textArgs
interp hide {} $w
# Install the alias:
interp alias {} $w {} graphareaCmd $w
foreach opt [array names options] {
$w configure $opt $options($opt)
}
return $w ;# the original canvas
}
proc graphareaCmd {self cmd args} {
switch -- $cmd {
configure {eval graphareaConfigure $self $cmd $args}
cget {eval graphareaCget $self $args}
showgrid { ;# draw a grid of lines on the graph
set x0 [$self cget xmin]
set x1 [$self cget xmax]
set y0 [$self cget ymin]
set y1 [$self cget ymax]
set dx [lindex $args 0]
set x [expr int($x0/$dx)*$dx] ;# choose lines s.t. zero is a line.
while {$x<$x1} {
$self addline $x $y0 $x $y1
set x [expr $x+$dx]
}
set dy [lindex $args 1]
set y [expr int($y0/$dy)*$dy]
while {$y<$y1} {
$self addline $x0 $y $x1 $y
set y [expr $y+$dy]
}
}
addline {
set args [eval concat $args]
set penc [$self cget pencolour];# get "local" pen colour name
set x0 [$self cget xmin]
set x1 [$self cget xmax]
set y0 [$self cget ymin]
set y1 [$self cget ymax]
set wid [$self cget -width]
set ht [$self cget -height]
foreach {x y} $args {
lappend xylist [expr $wid*double($x-$x0)/($x1-$x0)]
lappend xylist [expr $ht*double($y-$y1)/($y0-$y1)]
}
interp invokehidden {} $self create line $xylist -fill $penc
}
createcanvasoption { ;# add text (eg) positioned at scaled position
set args [eval concat $args] ;# this removes a level of curlies if necessary from the list.
set penc [$self cget pencolour];# get "local" pen colour name
set x0 [$self cget xmin]
set x1 [$self cget xmax]
set y0 [$self cget ymin]
set y1 [$self cget ymax]
set wid [$self cget -width]
set ht [$self cget -height]
lset args 1 [expr int($wid*double([lindex $args 1]-$x0)/($x1-$x0))]
lset args 2 [expr int($ht*double([lindex $args 2]-$y1)/($y0-$y1))]
switch [lindex $args 0] {
{text} {
interp invokehidden {} $self create text [lindex $args 1] [lindex $args 2] -fill $penc -text "[lindex $args 3]"
}
}
}
}
}
proc graphareaConfigure {self cmd args} {
# 3 scenarios:
#
# $args is empty -> return all options with their values
# $args is one element -> return current values
# $args is 2+ elements -> configure the options
switch [llength $args] {
0 { ;# return argument values
set result [ option get $self -pencolour ""]
# default options:
lappend result [interp invokehidden {} $self configure ]
# lappend result [uplevel 1 $self cconfigure]
return $result
}
1 {
switch -- $args {
{-pencolour} {return [uplevel 1 [list interp invokehidden {} $self configure -pencolour]]}
default {return [uplevel 1 interp invokehidden {} $self configure $args]}
}
}
default { ;# >1 arg - an option and its value
# go through each option:
foreach {option value} $args {
switch -- $option {
{-xmin} -
{-xmax} -
{-ymin} -
{-ymax} -
{-pencolour} {
global $self.[string range $option 1 end]
set $self.[string range $option 1 end] $value
}
default {puts " default $option, $value for $self";$self configure $option $value}
;#$self configure $option $value
}
}
return {}
}
}
}
proc graphareaCget {self args} {
# cget defaults done by the canvas cget command
#puts "In graphareaCget option $self $args"
if {[info exists ::$self.$args]} {upvar #0 $self.$args val; return $val}
switch -- $args {
{-pencolour} {return [uplevel 1 [list interp invokehidden {} $self cget pencolour]]}
default {return [uplevel 1 [list interp invokehidden {} $self cget $args]]}
}
}
if {1==0} {This section exercises 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 objecty is high. (EG Halley's comet spends 75 years orbiting the sun, of which less than 3 months is typically visible from earth.)
}
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}]]
}
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 ""
set pos4rk ""
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 $t [lindex $xvec 0] $t2 [lindex $result 0]
lappend pos4rk $t [lindex $xvec4 0] $t2 [lindex $res4 0]
set t $t2
set xvec $result
set xvec4 $res4
}
$dr2 configure -pencolour red
$dr2 addline $pos
interp invokehidden {} $dr2 create text 120 10 -fill red -text "Runge-Kutta-Fehlberg error correcting"
$dr2 configure -pencolour green
$dr2 addline $pos4rk
interp invokehidden {} $dr2 create text 120 24 -fill green -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 ""
while { $t < $tmax } {
set res4 [::math::calculus::rungeKuttaStep $t $tstep $xvec4 dampened_oscillator]
set t2 [expr {$t+$tstep}]
lappend pos $t [lindex $xvec4 0] $t2 [lindex $res4 0]
set t $t2
set xvec4 $res4
}
$dr2 configure -pencolour black
$dr2 addline $pos
interp invokehidden {} $dr2 create text 120 38 -fill black -text "Runge-Kutta short step"
set damp [expr {$damp+0.1}]
update
}An even better test is a particle bouncing between 2 walls - the forces are zero except when the particle hits the wall. The wall is stiff, so the force is very high briefly, requiring very short steps during the collision, but allowing long steps between collisions. The sample below correctly predicts the collisions at 1,3,5,7,9 and 11 seconds.
# wallhit force function: - 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 f 0
if {$x>1.0} { set f [expr -1.e7*($x-1.0)]
} elseif {$x<-1.0} { set f [expr -1.e7*($x+1.0)]
}
return [list $x1 $f]
}
set xvec { 0.0 1.0 } ;# initial position, velocity
set errvec {0.000001 1.e-8} ;# want accurate rebounds
set t 0.0
set tmax 12.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.2 ;# 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 addline $pos
See the reference [L1 ] for more details of RKF method.