>

prob1-4-4.mws

Although this worksheet is labeld as a solution to problem 1.4 #4, it can easily be modified to handle other problems in this section.

>

Gett ing started

Every Maple worksheet should begin by re-initializing the Maple "kernel" and loading the additional packages that we are most likely to use.

> restart;

> with( plots ):

> with( DEtools ):

>

Solution

The problem involves the following information

> MODEL := diff( y(t) , t ) = sin( y(t) ) ; # define the differential equation model

> IC := y(0) = 1; # specify the initial condition

>

Preliminary Analysis

We begin the analysis of this problem with the slope field.

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

> DOMAIN := t = 0 .. 3; # specify a reasonable interval for the independent variable

> RANGE := y = -0.5*Pi .. 1.5*Pi ; # specify a reasonable interval for the dependent variable

>

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

> plotSLOPE;

>

When the solution curve through the specified initial condition is added, the picture looks like

> plotSOLN:= DEplot( MODEL, VAR, DOMAIN, RANGE, [ [IC] ], linecolor = BLUE, arrows = MEDIUM ):

> plotSOLN;

Notice how the solution curve is always tangential to the slope field. (This is a good way to check that your "answer" is reasonable.)

>

Even though it is not asked for in this problem, let's see what Maple can do for an analytic solution.

> SOLN := dsolve( { MODEL, IC }, VAR );

UGH! No wonder we haven't spent more time worrying about explicit solutions!!!

>

Numerical Solution using Euler's Method

Now we get to the numerical solution. To do this we create a new command (or, actually, procedure), called Euler .

DO NOT CHANGE ANYTHING IN THE FOLLOWING INPUT REGION - YOU HAVE BEEN WARNED !!!

> Euler := proc( MODEL, IC, DOMAIN, N )

> local a, b, dt, dy, f, i, tt, yy, LISTpts, vars;

> vars := (op(2,lhs(MODEL) ), op([1,0],lhs(MODEL) ) );

> f := unapply( rhs(MODEL), vars );

> a := op(1,rhs(DOMAIN)); b := op(2,rhs(DOMAIN));

> dt := (b-a)/N;

> tt := a;

> yy := rhs(IC);

> LISTpts := [ tt, yy ];

> for i from 1 to N do

> dy := evalf( f(tt,yy) ) * dt;

> yy := yy + dy;

> tt := tt + dt;

> LISTpts := LISTpts, [ tt, yy ];

> od; RETURN( [ LISTpts ] );

> end:

>

Let's continue our preparations by checking that all necessary assignments have been made.

> MODEL;

> IC;

> VAR;

> DOMAIN:= t = 0 .. 8 ; # we've changed 3 to 8 for purposes of illustration -- play with it!

The last ingredient is the number of steps to be used in approximating the solution on the above domain. What value of delta t are we using here?

> Nsteps := 4;

The points, and corresponding plot, computed by Euler's method are:

> ptsE := Euler( MODEL, IC, DOMAIN, Nsteps );

> plotE := plot( ptsE, color = GREEN, thickness = 2 ):

> plotE;

>

Note that this plot uses a different vertical scale than the earlier plots. We could modify the plot command, but a better idea is to overlay the approximate solution on the slope field and see how well the solution follows the slope field. Since we have changed the domain we better recompute the slope field.

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

> display( [ plotSLOPE, plotE ] );

Observe that the approximate solution starts out following the slope field but starts to cross the slope field towards the end of each segment of the plotted solution.

>

We can also compare the Euler approximation with the "true solution" (which is itself, if truth be told, only a better approximation).

> plotSOLN:= DEplot( MODEL, VAR, DOMAIN, RANGE, [ [IC] ], linecolor = BLUE, arrows = MEDIUM ):

> display( [plotE, plotSOLN] );

>

Maple's built-in Euler method

The following much more concise code shows how to use Maple's built-in dsolver (with the classical forward Euler method) to get the same result. Unfortunately the results can be a little unpredictable.

> dt := (8-0)/Nsteps;

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

>

> NUMpts := 4;

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

>

>

> restart;