>

BDH2-5.mws

This worksheet is designed to support the reading of the text's discussion of the van der Pol equation (pp. 188 - 190). In particular, it contains a Maple implementation of Euler's method for 2-dimensional systems complete with examples of creating phase portraits and graphs of the (approximate) solutions.

>

Getting started

Defining the model

All differential equation models begin with a differential equation. If the problem is an intial value problem, an initial condition is also needed. Replace the question marks ( ? ) in the following input regions to define the relevant ODE model.

> ode1 := diff( x(t), t ) = y(t) ;

> ode2 := diff( y(t), t ) = (1-x(t)^2)*y(t) - x(t) ;

> MODEL := { ode1, ode2 } ; # define the differential equation model

> VAR := { x(t), y(t) } ; # identify the variables in the model

>

Direction Field, Phase Portrait, and Solution Curves ( DEplot )

Replace the question marks ( ? ) in the following input regions to create the direction field for your model.

> IC := [ x(0) = 1, y(0) = 1 ] ; # specify an initial condition (not needed for direction field)

> DOMAIN := t = 0 .. 2.5 ; # specify reasonable interval for indep var

> RANGE := x = -3 .. 3, y = -3 .. 3 ; # specify reasonable intervals for all dep vars

>

Direction Field

> plotDIR := DEplot( MODEL, VAR, DOMAIN, RANGE, arrows = MEDIUM,

> title = `Direction Field for van der Pol equation` ):

> plotDIR;

>

Phase Portrait

Remember that the stepsize option can be used to obtain a smoother curve if the default graph is too coarse.

> plotPHASE := DEplot( MODEL, VAR, DOMAIN, [ IC ], RANGE,

> scene=[ x, y ], arrows=NONE, linecolor=BLUE,

> title = `Phase Portrait for van der Pol equation` ):

> plotPHASE;

>

Solution Curves

> P1:= DEplot( MODEL, VAR, DOMAIN, [ IC ],

> scene=[ t, x ], arrows=NONE, linecolor=BLUE ):

> P2:= DEplot( MODEL, VAR, DOMAIN, [ IC ],

> scene=[ t, y ], arrows=NONE, linecolor=GREEN ):

> display( [ P1, P2 ] , title = `x and y vs. t` ); # combined solution curves

> #display( array([P1,P2]), title = `?` ); # side-by-side solution curves

>

Numerical Solutions ( dsolve,numeric )

We can also use this command to use a specified computational procedure. Most of these that are numerical (there are plenty of "symbolic" methods as well) require you to tell Maple how large the change in the independent variable should be (" stepsize "). So far the only method we know is Euler's, known to Maple as the "forward" Euler method from its " classical " collection. The main thing to investigate here is how the accuracy of the solution (and the time of computation) is affected by the stepsize. Let's recall what we are working with.

> MODEL;

> IC;

> VAR;

> DOMAIN;

>

Method # 1 - Euler2 (a custom-defined procedure for our uses)

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 := ?;

>

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 ) ]:

> plotXY := plot( ptsXY, color = GREEN, thickness = 2 ):

> plotXY;

>

Another way to view the approximate solution is to overlay the computed solution on the direction field (and Maple's phase portrait).

> display( [ plotDIR, plotPHASE, plotXY ] );

>

In a similar way, the two solution curves can be created:

> ptsTX := [ seq( [ pt[1], pt[2] ], pt=ptsE ) ]:

> plotTX := plot( ptsTX, color = GREEN, thickness = 2 ):

> plotTX;

>

> ptsTY := [ seq( [ pt[1], pt[3] ], pt=ptsE ) ]:

> plotTY := plot( ptsTY, color = GREEN, thickness = 2 ):

> plotTY;

>

And, again together with Maple's solution curve

> display( [ plotTX, P1 ] );

> display( [ plotTY, P2 ] );

>

Method #2 - dsolve,numeric and odeplot (built-in Maple commands) -- NOT YET UPDATED FOR SYSTEMS!

There are two steps in using Maple to work with numerically computed solutions to an initial value problem. First, the dsolve command is used with the numeric (and method ) option to obtain the solution.

>

Note that this command requires specification of the stepsize (not the number of steps). Fortunately, this conversion is not difficult to make.

> dt := ?

> solnE:= dsolve( { MODEL, IC }, VAR, numeric, method = classical[foreuler], stepsize = dt );

The second step is to decide how to use this solution. A common choice is to create a plot of the solution using odeplot (from the plots package). The last ingredient needed for this command is the number of points to use in the plot of the approximate solution.

> NUMpts := ?

> odeplot( solnE, [t, y(t)], rhs(DOMAIN), style = line, color = BLUE, numpoints = NUMpts );

>

Notes: