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
@time for n=1:npts-1
tn=(n-1)*h;
f1left=yH1[n]+h*f1(tn,yH1[n],yH2[n])
f2left=yH2[n]+h*f2(tn,yH1[n],yH2[n])
f1right=f1(tn+h,f1left,f2left)
f2right=f2(tn+h,f1left,f2left)
fH=(h/2)*(f1left+f1right)
yH1[n+1]=yH1[n]+fH
yH2[n+1]=yH2[n]+fH
end
Results can be shown in two ways:
plot([[y1[i] for i=1:npts], [y2[i] for i=1:npts]])
plot([(y1[i], y2[i]) for i=1:npts])
http://tutorial.math.lamar.edu/Classes/DE/Bernoulli.aspx
http://calculuslab.deltacollege.edu/ODE/7-C-3/7-C-3-h.html
Backward Euler method, or Implicit Euler method,