m |
m |
||
| Line 1: | Line 1: | ||
| + | We now turn to 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} | ||
| + | \frac{dx}{dt} &= \alpha x - \beta x y, \\ | ||
| + | \frac{dy}{dt} &= \delta x y - \gamma y, | ||
| + | \end{align} | ||
| − | + | The functions are implemented as: | |
<syntaxhighlight lang="python"> | <syntaxhighlight lang="python"> | ||
| − | + | function f1(t,y1, y2) | |
| − | + | α*y1-β*y1*y2 | |
| − | + | end | |
| + | function f2(t,y1, y2) | ||
| + | δ*y1*y2-γ*y2 | ||
end | end | ||
</syntaxhighlight> | </syntaxhighlight> | ||
| − | |||
Euler's method in this case reads: | Euler's method in this case reads: | ||
| Line 18: | Line 24: | ||
\end{align}$$ | \end{align}$$ | ||
| − | + | Heun's version reads: | |
<syntaxhighlight lang="python"> | <syntaxhighlight lang="python"> | ||
| − | + | @time for n=1:npts-1 | |
| − | + | y1[n+1]=y1[n]+(h/2)*(f1((n-1)*h,y1[n],y2[n])+f1(n*h,y1[n]+h*f1((n-1)*h,y1[n],y2[n]),y2[n])) | |
| − | + | y2[n+1]=y2[n]+(h/2)*(f2((n-1)*h,y1[n],y2[n])+f2(n*h,y1[n],y2[n]+h*f2((n-1)*h,y1[n],y2[n]))) | |
| − | + | ||
| − | + | ||
end | end | ||
</syntaxhighlight> | </syntaxhighlight> | ||
| + | |||
| + | RK4's version reads: | ||
| + | |||
| + | |||
<syntaxhighlight lang="python"> | <syntaxhighlight lang="python"> | ||
We now turn to 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
Euler's method in this case reads:
$$\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}$$
Heun's version reads:
@time for n=1:npts-1 y1[n+1]=y1[n]+(h/2)*(f1((n-1)*h,y1[n],y2[n])+f1(n*h,y1[n]+h*f1((n-1)*h,y1[n],y2[n]),y2[n])) y2[n+1]=y2[n]+(h/2)*(f2((n-1)*h,y1[n],y2[n])+f2(n*h,y1[n],y2[n]+h*f2((n-1)*h,y1[n],y2[n]))) end
RK4's version reads:
npts=10000 h=.001 y1=0.0*collect(1:npts); y2=0.0*collect(1:npts); y1[1]=1; y2[1]=.1; α=2/3; β=4/3; γ=1; δ=1; @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
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,