m |
m |
||
Line 103: | Line 103: | ||
$$ | $$ | ||
\begin{align} | \begin{align} | ||
− | k_1 &= h\vec f(t_n, \vec y_n)\,, \\ | + | t_{n+1} &= t_n + h\,, \\ |
− | k_2 &= h\vec f\left(t_n + \frac{h}{2}, \vec y_n + \frac{k_1}{2}\right)\,, \\ | + | \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)\,. |
− | k_3 &= h\vec f\left(t_n + \frac{h}{2}, \vec y_n + \frac{k_2}{2}\right)\,, \\ | + | \end{align} |
− | k_4 &= h\vec f\left(t_n + h, \vec y_n + k_3\right)\,, | + | $$ |
+ | |||
+ | 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{k_1}{2}\right)\,, \\ | ||
+ | \vec k_3 &= h\vec f\left(t_n + \frac{h}{2}, \vec y_n + \frac{k_2}{2}\right)\,, \\ | ||
+ | \vec k_4 &= h\vec f\left(t_n + h, \vec y_n + k_3\right)\,, | ||
\end{align} | \end{align} | ||
$$ | $$ |
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,