m |
m |
||
Line 1: | Line 1: | ||
− | We now turn to | + | 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 [https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations Volterra-Lotka model of prey-predator] interactions. This is modelled with the pair of coupled ODE: | ||
\begin{align} | \begin{align} | ||
Line 36: | Line 54: | ||
</syntaxhighlight> | </syntaxhighlight> | ||
+ | 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} | $$\begin{align} | ||
Line 54: | Line 74: | ||
Heun's version reads: | 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} | $$\begin{align} | ||
Line 60: | Line 84: | ||
\end{align}$$ | \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: | |
− | In code: | + | |
<syntaxhighlight lang="python"> | <syntaxhighlight lang="python"> | ||
@time for n=1:npts-1 | @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 | end | ||
</syntaxhighlight> | </syntaxhighlight> |
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
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=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
RK4's version reads:
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,