Crash Course in Scientific Computing
XVIII. Runge-Kutta
We now cover one of the most powerful predictor-corrector algorithms of all, one which is so accurate that most computer algorithms use it as default: the fourth-order Runge-Kutta Method. RK methods are optimized predictor-corrector methods, so they push the idea of Heun's method to the top. The RK4 method finds the optimum correction-prediction (in a way which we do not detail here) to involve these intermediate points:
\begin{align} k_1 &= h\ f(t_n, y_n)\,, \\ k_2 &= h\ f\left(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}\right)\,, \\ k_3 &= h\ f\left(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}\right)\,, \\ k_4 &= h\ f\left(t_n + h, y_n + k_3\right)\,, \end{align}
and then the calculation is explicit:
\begin{align} t_{n+1} &= t_n + h\,, \\ y_{n+1} &= y_n + \tfrac{1}{6}\left(k_1 + 2k_2 + 2k_3 + k_4 \right)\,. \end{align}
In code:
@time for n=1:npts-1 tn=(n-1)*h-t0; k1=h*f(tn,y[n]); k2=h*f(tn+(h/2),y[n]+k1/2); k3=h*f(tn+(h/2),y[n]+k2/2); k4=h*f(tn+h,y[n]+k3); y[n+1]=y[n]+(1/6)*(k1+2*k2+2*k3+k4) end
The error analysis on the same problems as above now yields a very different result:
herror=[]; for h in hlist npts=round(Int64,tmax/h); y=0.0*collect(1:npts); y[1]=-24/37; @time for n=1:npts-1 tn=(n-1)*h; k1=h*f(tn,y[n]); k2=h*f(tn+(h/2),y[n]+k1/2); k3=h*f(tn+(h/2),y[n]+k2/2); k4=h*f(tn+h,y[n]+k3); y[n+1]=y[n]+(1/6)*(k1+2*k2+2*k3+k4) end append!(herror,mean(abs.([-(24/37)*cos(3t)-(4/37)*sin(3t) for t in ((1:npts).-1).*h]-y))) end
0.006093 seconds (3.91 k allocations: 64.688 KiB) 0.007069 seconds (41.92 k allocations: 674.766 KiB) 0.071775 seconds (437.74 k allocations: 6.871 MiB) 0.805397 seconds (4.40 M allocations: 68.998 MiB, 1.76% gc time) 6.734973 seconds (43.98 M allocations: 690.260 MiB, 1.41% gc time) 77.495034 seconds (439.82 M allocations: 6.741 GiB, 1.57% gc time)
scatter(herror, xlabel="time step", ylabel="mean abs error", xticks=(1:6, string.(hlist)), legend=false, yscale=:log10) plot!(x->10^4*herror[1]*10.0^-(4x), 1, 6, ls=:dot, lw=3, ylims=(10^-15,10^-3))
The gain from $h=0.001$ to $h=0.0001$ is marginally better, and for smaller still $h$, it actually becomes worse, due to rounding errors. It works so well that many, if not most, resources would use the RK4 to integrate ODE.
We now turn to the case of coupled equations, in which case the method reads
$$\vec y'=\vec f(t,\vec y)$$
where
$$ \vec y\equiv \begin{pmatrix} y_1\\y_2 \end{pmatrix} \quad\text{and}\quad \vec f\equiv \begin{pmatrix} f_1\\f_2 \end{pmatrix}\,. $$
The methods are the same as before just applying directly to the vectors instead. An interesting coupled system of ODE, namely, the Volterra-Lotka model of prey-predator interactions. This is modelled with the pair of coupled ODE:
\begin{align} \frac{dx}{dt} &= \alpha x - \beta x y, \\ \frac{dy}{dt} &= \delta x y - \gamma y, \end{align}
The functions are implemented as:
function f1(t,y1, y2) α*y1-β*y1*y2 end function f2(t,y1, y2) δ*y1*y2-γ*y2 end
Some common initialization for all methods, in particular, we now need two arrays of values to store the evolution of $x$ (in y1) and $y$ (in y2):
tmax=100 h=.001 npts=round(Int,tmax/h) y1=0.0*collect(1:npts); y2=0.0*collect(1:npts); y1[1]=1; y2[1]=.1; α=2/3; β=4/3; γ=1; δ=1;
Euler's method in this case reads:
$$\vec y_{n+1}=\vec y_n+h\vec f(t_n,\vec y_n)$$
or, breaking down this equation componentwise:
$$\begin{align} y_{1,n+1}&=y_{1,n}+hf_1(t_n,y_{1,n},y_{2,n})\\ y_{2,n+1}&=y_{2,n}+hf_2(t_n,y_{1,n},y_{2,n}) \end{align}$$
In code:
@time for i=1:npts-1 y1[i+1]=y1[i]+h*f1((i-1)*h,y1[i],y2[i]) y2[i+1]=y2[i]+h*f2((i-1)*h,y1[i],y2[i]) end
Heun's version reads:
$$\vec y_{n+1}=\vec y_n+{h\over 2}\left(\vec f(t_n,\vec y_n)+\vec f(t_{n+1},\vec y_n+h\vec f(t_i,\vec y_n)\right)$$
or, breaking down this equation componentwise:
$$\begin{align} y_{1,n+1} &= y_{1,n} + {h\over2} [f(t_n, y_{1,n}, y_{2,n}) + f(t_n + h, y_{1,n} + h f(t_n, y_{1,n}, y_{2,n}), y_{2,n})]\\ y_{2,n+1} &= y_{2,n} + {h\over2} [f(t_n, y_{1,n}, y_{2,n}) + f(t_n + h, y_{1,n}, y_{2,n} + h f(t_n, y_{1,n}, y_{2,n}))] \end{align}$$
In code; where this time we'll find both more convenient and efficient to define auxiliary quantities, that are furthermore invoked more than once:
@time for n=1:npts-1 tn=(n-1)*h; f1left=f1(tn,yH1[n],yH2[n]) f2left=f2(tn,yH1[n],yH2[n]) f1right=f1(tn+h,yH1[n]+h*f1left,yH2[n]+h*f2left) f2right=f2(tn+h,yH1[n]+h*f1left,yH2[n]+h*f2left) yH1[n+1]=yH1[n]+(h/2)*(f1left+f1right) yH2[n+1]=yH2[n]+(h/2)*(f2left+f2right) end
RK4's version reads:
$$ \begin{align} t_{n+1} &= t_n + h\,, \\ \vec y_{n+1} &= \vec y_n + \tfrac{1}{6}\left(\vec k_1 + 2\vec k_2 + 2\vec k_3 + \vec k_4 \right)\,. \end{align} $$
in terms of the intermediate quantities:
$$ \begin{align} \vec k_1 &= h\vec f(t_n, \vec y_n)\,, \\ \vec k_2 &= h\vec f\left(t_n + \frac{h}{2}, \vec y_n + \frac{\vec k_1}{2}\right)\,, \\ \vec k_3 &= h\vec f\left(t_n + \frac{h}{2}, \vec y_n + \frac{\vec k_2}{2}\right)\,, \\ \vec k_4 &= h\vec f\left(t_n + h, \vec y_n + \vec k_3\right)\,, \end{align} $$
In code
@time for n=1:npts-1 tn=(n-1)*h; # Intermediate steps k11=h*f1(tn,yRK1[n],yRK2[n]); k12=h*f2(tn,yRK1[n],yRK2[n]); k21=h*f1(tn+(h/2),yRK1[n]+k11/2,yRK2[n]+k12/2); k22=h*f2(tn+(h/2),yRK1[n]+k11/2,yRK2[n]+k12/2); k31=h*f1(tn+(h/2),yRK1[n]+k21/2,yRK2[n]+k22/2); k32=h*f2(tn+(h/2),yRK1[n]+k21/2,yRK2[n]+k22/2); k41=h*f1(tn+h,yRK1[n]+k31,yRK2[n]+k32); k42=h*f2(tn+h,yRK1[n]+k31,yRK2[n]+k32); # Real step yRK1[n+1]=yRK1[n]+(1/6)*(k11+2*k21+2*k31+k41); yRK2[n+1]=yRK2[n]+(1/6)*(k12+2*k22+2*k32+k42); end
Results can be shown in two ways:
1 - As the respective solutions. This shows the oscillations of predators and preys
pRK=plot([[yRK1[i] for i=1:240], [yRK2[i] for i=1:240]], lw=3, xlabel="time", ylabel="populations", label="RK4")
2 - As a phase-space solution, showing the simultaneous state of the system (here, ecosystem), with time as a parametric parameter:
psRK=plot([(yRK1[i], yRK2[i]) for i=1:npts], lw=1, xlabel="x", ylabel="y", legend=false)
This comparing the three methods together (the various plots have been evaluated as pE, pH, etc.):
plot(pE,pH,pRK,psE,psH,psRK,layout=(2,3),dpi=150)
The time-evolutions at the top are over considerably fewer cycles for Heun & RK4 methods than those simulated to produce the phase space solutions, bottom. You can see how the RK4 method stays stable (which is the correct behaviour) while the Heun method slowly drifts towards the stable point.