Euler's method for 2D systems

October 7, 1998

Here is how use Maple to perform Euler's method for a 2D system.
It's a straightfoward extension of the 1D version.
We'll apply it to the Van der Pol initial value problem:

>

dt:=0.50: n:=ceil(4.0/dt):   t:=0.0:   x:=1.0:   y:=1.0:   txdata[0] := [t,x]:   tydata[0] := [t,y]:   xydata[0] := [x,y]: txydata[0] := [t,x,y]:   for i from 1 to n do        mk := y:        nk := -x + (1-x^2)*y:        x := x + dt*mk:        y := y + dt*nk:        t := t + dt:        txdata[i] := [t,x]:      tydata[i] := [t,y]:     xydata[i] := [x,y]:     txydata[i] := [t,x,y]: od:   t:='t': x:='x': y:='y':                      # (clean-up) txdatalist := [ seq( txdata[i], i=0..n ) ]:   tydatalist := [ seq( tydata[i], i=0..n ) ]:   xydatalist := [ seq( xydata[i], i=0..n ) ]: txydatalist := [ seq( txydata[i], i=0..n ) ]: txymatrix := convert(txydatalist,matrix);

Here is a way to plot the and functions.

>	plot( {txdatalist,tydatalist}, `t`=0..4, `x and y`=-3.5..1.8);

Here is how to get a phase plane plot.

>	plot( xydatalist );

>	with(DEtools): with(plots):

Here the phase plane plot is combined with the slope field plot.

>	fieldplot := dfieldplot([D(x)(t)=y(t),D(y)(t)=-x(t)+(1-x(t)^2)*y(t)], [x(t),y(t)],t=0..10,x=-3..3,y=-5..5): eulerplot := plot( xydatalist,color=blue): display([fieldplot,eulerplot]);

>

dt:=0.05: n:=ceil(4.0/dt):   t:=0.0:   x:=1.0:   y:=1.0:   txdata[0] := [t,x]:   tydata[0] := [t,y]:   xydata[0] := [x,y]: txydata[0] := [t,x,y]:   for i from 1 to n do        mk := y:        nk := -x + (1-x^2)*y:        x := x + dt*mk:        y := y + dt*nk:        t := t + dt:        txdata[i] := [t,x]:      tydata[i] := [t,y]:     xydata[i] := [x,y]:     txydata[i] := [t,x,y]: od:   t:='t': x:='x': y:='y':                      # (clean-up) txdatalist2 := [ seq( txdata[i], i=0..n ) ]:   tydatalist2 := [ seq( tydata[i], i=0..n ) ]:   xydatalist2 := [ seq( xydata[i], i=0..n ) ]: txydatalist2 := [ seq( txydata[i], i=0..n ) ]: txymatrix2 := convert(txydatalist2,matrix);

>	eulerplot2 := plot( xydatalist2,color=green): display([fieldplot,eulerplot,eulerplot2]);

>

dt:=0.005: n:=ceil(4.0/dt):   t:=0.0:   x:=1.0:   y:=1.0:   txdata[0] := [t,x]:   tydata[0] := [t,y]:   xydata[0] := [x,y]: txydata[0] := [t,x,y]:   for i from 1 to n do        mk := y:        nk := -x + (1-x^2)*y:        x := x + dt*mk:        y := y + dt*nk:        t := t + dt:        txdata[i] := [t,x]:      tydata[i] := [t,y]:     xydata[i] := [x,y]:     txydata[i] := [t,x,y]: od:   t:='t': x:='x': y:='y':                      # (clean-up) txdatalist3 := [ seq( txdata[i], i=0..n ) ]:   tydatalist3 := [ seq( tydata[i], i=0..n ) ]:   xydatalist3 := [ seq( xydata[i], i=0..n ) ]: txydatalist3 := [ seq( txydata[i], i=0..n ) ]: txymatrix3 := convert(txydatalist3,matrix): eulerplot3 := plot( xydatalist3,color=black): display([fieldplot,eulerplot,eulerplot2,eulerplot3]);

Below we list the endpoints of the second ( ) and third ( ) approximations and the difference of the two.

>	txydatalist3[nops(txydatalist3)];

>	txydatalist2[nops(txydatalist2)];

>	txydatalist3[nops(txydatalist3)]-txydatalist2[nops(txydatalist2)];

Now, Maple has even more powerful numerics, using not Euler's Method but a Runge-Kutta approximation of degree 4.

>	g:=dsolve({D(x)(t)=y(t),D(y)(t)=-x(t)+(1-x(t)^2)*y(t),x(0)=1,y(0)=1}, [x(t),y(t)],type=numeric);

Here's what this gets for an endpoint.

>

g(4);

>	with(plots):

And here's the and graphs.

>	plot1:=odeplot(g,[t,x(t)],0..4): plot2:=odeplot(g,[t,y(t)],0..4,color=blue): display([plot1,plot2]);

And the phase plane.

>	odeplot(g,[x(t),y(t)],0..4);

And here we have them all together, the latest addition in cyan, barely visible next to the black ( ) curve.

>	dsolveplot:=odeplot(g,[x(t),y(t)],0..4,color=cyan): display([fieldplot,eulerplot,eulerplot2,eulerplot3,dsolveplot]);

Here's the difference between the dsolve,numeric endpoint and the Euler's Method ( ) endpoint.

>	subs(g(4),[t,x(t),y(t)])-txydatalist3[nops(txydatalist3)];

These two are fairly close.  No wonder the cyan curve is barely visible.