Implementation of Euler's Method for systems

Rabbits and Wolves

ODE:  

r’  =  3r - wr/2
w’ = -w/10 + 7rw/10

vf[r_,w_]={3 r-w r/2,-w/10+7 r w/10}

df=VectorPlot[vf[r,w]/Norm[vf[r,w]],{r,0,4},{w,0,10},VectorStyle->Arrowheads[.03]]

Equilibria: (0,0), (1/7,6)

Show[df,ListPlot[{{0,0},{1/7,6}},PlotStyle->Red]]

dfcloseup1=VectorPlot[vf[r,w]/Norm[vf[r,w]],{r,-0.004,0.004},{w,-0.004,0.004},VectorStyle->Arrowheads[.03],VectorPoints->Fine]

VectorPlot[vf[r,w]/Norm[vf[r,w]],{r,-0.04,0.04},{w,-0.004,0.004},VectorStyle->Arrowheads[.03],VectorPoints->Fine]

Show[dfcloseup2=VectorPlot[vf[r,w]/Norm[vf[r,w]],{r,1/7-0.01,1/7+0.01},{w,6-0.01,6+0.01},VectorStyle->Arrowheads[.03]],ListPlot[{{1/7,6}},PlotStyle->Red]]

Well, both of these are quite unsatisfying pictures.  I’ll guess the nontrivial equilibrium point (1/7,6) is a spiral sink.  Let’s try to see if we can confirm that with a solution.

soln[rzero_?NumericQ,wzero_?NumericQ]:=NDSolve[{dr'[t]==3 dr[t]-0.5 dr[t] dw[t],dw'[t]==-0.1 dw[t]+0.7 dw[t] dr[t],dr[0]==rzero,dw[0]==wzero},{dr,dw},{t,0,100}]
rr[t_?NumericQ,rzero_?NumericQ,wzero_?NumericQ]:=dr[t]/.soln[rzero,wzero][[1]]
ww[t_?NumericQ,rzero_?NumericQ,wzero_?NumericQ]:=dw[t]/.soln[rzero,wzero][[1]]

orbit1=ParametricPlot[{rr[s,1,1],ww[s,1,1]},{s,0,100},PlotStyle->{Thick,Red}];

Show[df,orbit1]

orbit2=ParametricPlot[{rr[s,1/7 +0.1,6],ww[s,1/7+0.1,6]},{s,0,100},PlotStyle->{Thick,Red}];

Show[df,orbit2]

My plots are inconclusive.  One can see what turns up via linearization.  If there really are periodic orbits, then I’ll guess we have a Hamiltonian.  In any case, let’s go ahead and try to implement Euler’s method.  Since the ODE is autonomous, I’ll take tzero = 0 and start with an arbitrary initial population distribution.  (Why does this make sense?)  I’ll also restrict to the time interval [0,10].  Finally, I’ll suppress the time, and compute the image in phase space.

euler[xzero_?NumericQ,yzero_?NumericQ,N_?NumericQ]:= For[j=1; image={{xzero,yzero}},j< N+1, j++,image=Append[image,image[[j]]+vf[image[[j]][[1]],image[[j]][[2]]] (10/N)]]

euler[1,1,2]

image

show2=Show[df,ListPlot[image,PlotStyle->Red],PlotRange->All]

euler[1,1,3]

image

show3=Show[df,ListPlot[image,PlotStyle->Red],PlotRange->{{0,15},{0,10}}]

euler[1,1,4]

image

show4=Show[df,ListPlot[image,PlotStyle->Red],PlotRange->{{0,15},{0,10}}]

euler[1,1,5]

image

show5=Show[df,ListPlot[image,PlotStyle->Red],PlotRange->{{0,15},{0,10}}]

euler[1,1,20]

showend=Show[df,ListPlot[image,PlotStyle->Red],PlotRange->{{0,15},{0,10}}]