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.
gwm 17-05-05 correction to recommended step length improves speed, and accuracy.
# 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
# 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} {#puts "Miniisiung $spp, expr $spp*$tstep"
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
# puts "[info level] at $t $tstep rec $stepest dt is [expr {double($tstep)/int(ceil( $tstep/$stepest))}] [lindex $result 0] [lindex $result 1]"
set i 0
set tsu $t
while {$i<$nsubsteps} {
set result [rungeKuttaFehlbergStep $tsu $dtsub $result $func errv ]
# puts "[info level] at $t dt is $dtsub [lindex $result 0] [lindex $result 1]"
# 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} {This section is the graph drawing widget option found under Another Graph 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
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 create line $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 create line $x0 $y $x1 $y
set y [expr {$y+$dy}]
}
}
create { ;# replaces the create default of a canvas which works in pixels
#- adds text, rectangle, oval, polygon... (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
# scale factor for draw space to pixels
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]
set idx 1 ;# where to start in args for coordinates
while {$idx<[llength $args]} {
if {[string is double [lindex $args $idx]]} {
lappend xylist [expr {int($wid*double([lindex $args $idx]-$x0)/($x1-$x0))}]
set args [lreplace $args $idx $idx]
lappend xylist [expr {int($ht*double([lindex $args $idx]-$y1)/($y0-$y1))}]
set args [lreplace $args $idx $idx]
} else {
incr idx
}
}
switch [lindex $args 0] {
{line} -
{text} {
lappend command [lindex $args 0] $xylist -fill $penc
if {[llength $args]>1} { set command [concat $command [lrange $args 1 end]] }
eval interp invokehidden {{}} $self create $command
}
{default} {
lappend command [lindex $args 0] $xylist -outline $penc
if {[llength $args]>1} { set command [concat $command [lrange $args 1 end]] }
eval interp invokehidden {{}} $self create $command
}
}
}
}
}
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} {
if {1==0} {This is the part I dont like - I create a global variable for each added option.
}
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 {
{-xmin} -
{-xmax} -
{-ymin} -
{-ymax} -
{-pencolour} {puts "Cget option $args"
return [uplevel 1 [list interp invokehidden {} $self cget [string range $args 1 end]]]}
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 object is high. (EG Halley's comet spends 76 years orbiting the sun, for less than 3 months is it 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 ""
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
}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.
puts "Type is $type"
# 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 1.0 } ;# 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.05 ;# 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 $posSee the reference [L1 ] for more details of RKF method.