Runge-Kutta methods of solving ODE Ordinary Differential Equations are numerical techniques to find the solution to an ODE. R-K methods are more stable and accurate than Euler Methods
A number of steps are used to advance a parameter through time (or space, depending on the equation). For example:
d/dt(N)=-k.N
(rate of change of number proportional to the population).
Euler method would try:
loop (until fed up) {
N=N-k.N*dt
time=time+dt
}Which deviates from the true solution very quickly since the population is changing during the time step; very short time steps will reduce the error but take a long time.
R-K methods use a numerical integration procedure which takes several (often 4) intermediate calculations to find an approximate integral of the equation. See reference 1 below for more details. The following code will solve d(N)/dt=-2*n with good accuracy.
# tcl runge-kutta solution of ODES # single variable - use rk4step (included for simplicity)
proc decay { time xn } { ;# simple decay proportional to amount (gives radioactive half life, eg)
return [expr -2.0*$xn]
}
proc rk4step {x time dt dxdtproc } { ;# proc is such that dx/dt = dxdtproc(t)
set k1 [$dxdtproc $time $x] ;# k1 = f (tn, xn)
set k2 [$dxdtproc [expr $time+($dt*0.5)] [expr $x + $dt*0.5*$k1]] ;# k2 = f (tn + h?2, xn + h?2 k1)
set k3 [$dxdtproc [expr $time +$dt*0.5] [expr $x + $dt*0.5*$k2] ] ;# k3= f (tn + h?2, xn + h?2 k2)
set k4 [$dxdtproc [expr $time +$dt] [expr $x + $dt*$k3] ] ;# k4= f (tn + h, xn + h k3)
return [expr $x + $dt/6.0 *($k1 + 2*$k2 + 2*$k3 + $k4)] ;# h/6(simpsons rule)
}
proc test {} {
set y 0
set mass 1
set dt 0.1 ;# try shortening dt to see numerical inaccuracy
# use single equation R_K scheme - log decay of radioactivity eg (dN/dt=-n-->> N(t)=exp(-t)
for {set time 0} {[expr $time<=2] } { set time [expr $time+$dt]} {
puts "At time $time $mass"
set mass [rk4step $mass $time $dt decay] ;# dmass/dt = -mass - mass=m0(1-exp(time/k))
}
}
testat time =0.5 decay should be exp(-1)=0.367879441171 this gives At time 0.5 0.367885238126 at time =1 decay should be exp(-2)=0.135335283237 this gives At time 1.0 0.135339548431
This reference is simple: [L1 ] (and useful).
The pdf reference is very comprehensive [L2 ] including adaptive stepsize control and higher order methods for error limitation.