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 solverred 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 [1] for more details of RKF method.[Category Mathematics|Category Physics|Category Science|Category Numerical Analysis]