Euler2 (a custom-defined procedure for two-dimensional systems)
DO NOT CHANGE ANYTHING IN THE FOLLOWING INPUT REGION - YOU HAVE BEEN WARNED !!!
> Euler2 := proc( MODEL, IC, DOMAIN, N )
> local a, b, dt, dx, dy, f, g, i, ode1, ode2, tt, xx, yy, LISTpts, vars;
> ode1 := op(1,MODEL);
> ode2 := op(2,MODEL);
> vars := (op(2,lhs(ode1) ), op([1,0],lhs(ode1) ), op([1,0],lhs(ode2) ) );
> f := unapply( rhs(ode1), vars );
> g := unapply( rhs(ode2), vars );
> a := op(1,rhs(DOMAIN)); b := op(2,rhs(DOMAIN));
> dt := (b-a)/N;
> tt := a;
> xx := subs( IC, vars[2](a) );
> yy := subs( IC, vars[3](a) );
> #yy := rhs(IC);
> LISTpts := [ tt, xx, yy ];
> for i from 1 to N do
> dx := evalf( f(tt,xx,yy) ) * dt;
> dy := evalf( g(tt,xx,yy) ) * dt;
> xx := xx + dx;
> yy := yy + dy;
> tt := tt + dt;
> LISTpts := LISTpts, [ tt, xx, yy ];
> od;
> RETURN( [ LISTpts ] );
> end:
DO NOT CHANGE ANYTHING IN THE PRECEDING INPUT REGION - YOU HAVE BEEN WARNED !!!
>
Let's recall what we are working with.
> MODEL;
> IC;
> VAR;
> DOMAIN ;
The number of steps to be used in Euler's method is
> Nsteps := 6;
>
The points generated by Euler's method are
> ptsE := Euler2( MODEL, IC, DOMAIN, Nsteps );
>
These data points can be used to create an approximate phase portrait as follows:
> ptsXY := [ seq( [ pt[2], pt[3] ], pt=ptsE ) ]:
> XYplot := plot( ptsXY, color = GREEN, thickness = 2 ):
> XYplot;
>
Another way to view the approximate solution is to overlay the computed solution on the direction field (and Maple's phase portrait).
> display( [ DIRplot, PHASEplot, XYplot ] );
>
Now, increase the time domain to show the solution for t=15. Recompute the Euler solution with 6 steps (dt=2.5). Next, compute the Euler solution with dt=0.25; how do these two solutions compare?