Julia programming

Revision as of 15:42, 17 March 2020

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

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:

$$\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:

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,