**Maple Tools for Differential Equations**

**A. J. Meir**

Copyright (C) A. J. Meir. All rights reserved.

This worksheet is for educational use only. No part of this publication may be reproduced or transmitted for profit in any form or by any means, electronic or mechanical, including photocopy, recording, or any information storage and retrieval system without prior written permission from the author. Not for profit distribution of the software is allowed without prior written permission, providing that the worksheet is not modified in any way and full credit to the author is acknowledged.

(This document is largely based on a similar document originally written by P. G. Schmidt.)

Maple will find
*symbolic solutions*
(that is,
*closed-form solutions*
or
*solution formulas*
) for many types of d.e.s and i.v.p.s, including separable and linear first-order equations. Moreover, Maple can compute and graph
*numerical solutions*
(that is,
*approximate solutions*
) for virtually all i.v.p.s of practical interest.

**1. Finding and Plotting Symbolic Solutions.**
To find a general solution for a given d.e., use Maple's
dsolve
command, whose general format is

dsolve(DE,unknown_function) .

The following command, for example, finds a general solution for the d.e. ' = (the radioactive-decay equation), where the unknown function is and is a parameter:

`> `
**dsolve(D(x)(t)=-k*x(t),x(t));**

Note that inside the dsolve command the derivative, ', should be written as and that the independent variable, , must be carried along (so that the d.e. becomes instead of ' = ). In Maple's response, denotes an arbitrary constant. (In general, Maple renders arbitrary constants as , , , and so forth.) In standard calculus notation, we would write the solution formula as

.

The next command produces a general solution for the d.e. '' + (the elastic-spring equation with ), where the unknown function is again and is a parameter:

`> `
**dsolve(D(D(x))(t)+a^2*x(t),x(t));**

Note that the second derivative, '', should be written as and that the independent variable must be carried along, just as in the previous example.

With a minor modification, we can use the dsolve command to find symbolic solutions of i.v.p.s; the modified format is

dsolve({DE,initial_condition(s)},unknown_function) .

Note the "curly braces" around the d.e. and initial condition(s). For example, let's solve the first-order i.v.p. ' = , :

`> `
**dsolve({D(x)(t)=-k*x(t),x(0)=x0},x(t));**

Next, we solve the second-order i.v.p. '' + , , x'(0) = :

`> `
**dsolve({D(D(x))(t)+a^2*x(t)=0,x(0)=x0,D(x)(0)=v0},x(t));**

In order to
*plot*
a solution generated by the
dsolve
command, we first have to give it a name, say,
. The following command sets
equal to the right-hand side of the immediately preceding solution formula:

`> `
**s:=rhs(%);**

Note: The percent (%) is just a placeholder for the
*last output generated before execution of the present command*
(which, in this case, happens to be the solution formula produced by the preceding
dsolve
command). Instead of using the percent, we could have given a name (say,
) to the solution formula and used that name in the
rhs
command:

`> `
**sol:=dsolve({D(D(x))(t)+a^2*x(t)=0,x(0)=x0,D(x)(0)=v0},x(t));**

`> `
**s:=rhs(sol);**

If we now specify the value of the parameter and the initial values and , we can plot just like any other function. Set, for example, , , and :

`> `
**a:=1; x0:=2; v0:=-3;**

In this case, the solution of the i.v.p. is

`> `
**s;**

and a graph can be generated with the usual plot command:

`> `
**plot(s,t=-10..20);**

To obtain the solution and its graph for other values of the parameter or other initial values and , just re-execute the last three commands, changing only the first one as needed. For further details on plotting, especially on simultaneous plotting of several graphs, see the handout " More on Basic Maple Commands".

**2. Plotting Direction Fields and Numerical Solutions.**
The basic Maple command for the graphical representation of direction fields and approximate solution curves is
DEplot
; it is part of a package of routines called
DEtools
, which must be explicitly loaded with the command

`> `
**with(DEtools):**

To plot a
*direction field*
for a first-order d.e. of the form
' =
, in the rectangle
,
, use the
DEplot
command in the form

DEplot(D(x)(t)=f(t,x(t)),x(t),t=a..b,x=c..d) .

The following command, for example, generates a direction field for ' = in the rectangle , :

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-6..6);**

The default number of grid points used is 20 in either coordinate direction. You can change that, if needed, by adding the option

dirgrid=[m,n]

inside the DEplot command, where m and n denote the desired number of grid points in the horizontal and vertical direction, respectively. Example:

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-6..6,dirgrid=[15,30]);**

In addition to a direction field, Maple will compute and display one or more
*approximate solution curves*
of a given d.e. if suitable initial conditions are specified in the form

{[x(t0)=x0],[x(t1)=x1],[x(t2)=x2],...} ,

immediately following the range specifications in the DEplot command. Note that the set of initial conditions is enclosed in curly braces, while each individual condition appears in brackets. Both, brackets and braces, must be present even if only a single initial condition is given.

For example, the following command produces a direction field for ' = , as before, along with the (approximate) graph of a particular solution, satisfying the initial condition :

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-6..6,{[x(0)=1]});**

Now we will add a second solution curve, satisfying the initial condition :

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-6..6,{[x(0)=1],[x(0)=-1]});**

In general, the given initial points do not have to lie inside the plot window.

If initial conditions are given, specification of the vertical range ( -range) is optional; if omitted, DEplot will choose the -range automatically, depending on where the solution curves "go". This is sometimes useful, but not always. Let's try it with the above example (same command as before, but without -range specification):

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-6..6,{[x(0)=1],[x(0)=-1]});**

Can you explain the poor outcome? (Hint: The
*exact*
solution with initial condition
has a vertical asymptote at
.)

If only the solution curves are of interest, display of the direction field arrows can be suppressed by adding the option

arrows=none

inside the DEplot command, following the specification of initial conditions. Example:

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-6..6,{[x(0)=1],[x(0)=-1]},arrows=none);**

Many further options are available. In fact, most of the options available with the basic plot command will also work with DEplot . Ver y useful for printing purposes is the option

linecolor=black ,

which causes the solution curves to be drawn in black rather than yellow (the default, which looks pitiful on a black-and-white printout). The color and appearance of the direction field arrows can be customized, too, but this will hardly ever be necessary. For details, consult the relevant help pages (by typing ?DEplot or ?plot,options).

Note: Options may be added to the DEplot command in any order, but must be preceded by the range specification(s) and initial condition(s) (if any). Example:

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-6..6,{[x(0)=1],[x(0)=-1]},
dirgrid=[10,10],linecolor=black,color=green,title=`dx/dt=x^2`);**

The default method that DEplot uses to compute approximate solutions of the given d.e. is "rk4", a so-called 4th order Runge-Kutta method, much more sophisticated than the classical Euler method. You can force Maple to use the Euler method by specifying the option

method=classical .

Both "rk4" and "Euler" employ a fixed stepsize; the default stepsize is ( )/20, where is the horizontal range ( -range) of the plot. Frequently, it will be necessary to reduce the stepsize, for improved numerical accuracy. This can be done by specifying the option

stepsize=h ,

where h denotes the desired stepsize, smaller than ( )/20 . Let's try the following example, where we approximate the same solutions of ' = as before, but with a larger vertical range ( ) and using the classical Euler method (instead of "rk4"):

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-9..9,{[x(0)=1],[x(0)=-1]},method=classical,linecolor=black);**

The horizontal part of the lower solution curve is clearly erroneous -- it certainly doesn't match the direction field! The (default) stepsize used in the plot is .3 (namely, ( )/20 with and ), and apparently this is not accurate enough. Reducing the stepsize should remedy the problem:

`> `
**DEplot(D(x)(t)=x(t)^2,x(t),t=-3..3,x=-9..9,{[x(0)=1],[x(0)=-1]},method=classical,stepsize=.1,linecolor=black);**

To accurately track a rapidly growing or highly oscillatory solution, it may be necessary to employ a very small stepsize and compute millions of data points. Not only will this consume a lot of time -- storing millions of data points may also tax Maple's memory. Fortunately, it is not nessary to
*store*
all the computed data points (the number of points used in creating the plot is anyway limited by the resolution of computer screen and laser printer). With the option

iterations=n ,

you can instruct Maple to store and plot only every n th of the computed data points. This option, if applied judiciously, can save a lot of time and memory, with no loss of numerical accuracy. As a drastic example, we will compute and display an approximate solution of the i.v.p.

, ,

on the interval
, using "rk4" (the default method) with a stepsize of .001 and 100 iterations between plotted points. (Thus, Maple will
*compute*
15,000 data points, but
*store and plot*
only 150.)

`> `
**DEplot(D(x)(t)=1-exp(t)*x(t)^2,x(t),t=0..15,{[x(0)=0]},
stepsize=.001,iterations=100,arrows=none,linecolor=black);**

Generating this plot still takes quite some time, but without the option iterations=100 , it would take about ten times as long, without any visible improvement! On the other hand, if you tried to speed things up by using a stepsize significantly greater than .001, the plot would not be reliable. Check it out!

Beyond "rk4" (the default method) and "Euler" ( method=classical ), you have a choice of several, much more sophisticated d.e. solvers. One of them is "rkf45", a so-called 4th/5th order Runge-Kutta-Fehlberg method, which can be invoked by specifying

method=rkf45

in the
DEplot
command. While "rk4" and "Euler" use a
*fixed*
stepsize, "rkf45" employs a
*variable*
stepsize that's automatically adjusted in accordance with the "stiffness" of the problem. As an adaptive (self-correcting) method, "rkf45" is more reliable than "rk4" or "Euler", but usually slower. But in "stiff" problems (like the previous example), "rkf45" offers a double advantage: You don't have to worry about the stepsize,
*and*
it is faster than "rk4" or "Euler"! (In general, controlling and adjusting the stepsize as needed is more efficient than using a very small stepsize throughout.)

Since "rkf45" controls the (numerical) stepsize automatically, specifying a stepsize when using this method has no effect on the numerical accuracy of the procedure; it does, however, affect the number of points plotted and, thereby, the apparent smoothness of the graph. More precisely, in conjunction with "rkf45", the specified stepsize becomes simply the horizontal distance between consecutive points plotted. For example, if we solve the previous problem with "rkf45", we should specify a stepsize of .1 to obtain the same number of points plotted as before, when we used "rk4" with a (numerical) stepsize of .001 and 100 iterations between points plotted. Let's try it:

`> `
**DEplot(D(x)(t)=1-exp(t)*x(t)^2,x(t),t=0..15,{[x(0)=0]},
method=rkf45,stepsize=.1,arrows=none,linecolor=black);**

DEplot
will not only solve
*first-order*
, but also
*higher-order*
i.v.p.s (although no direction fields can be generated for the latter). The following command, for example, produces an approximate solution of the i.v.p.
'' +
,
, x'(0) = -3, whose symbolic solution is
:

`> `
**DEplot(D(D(x))(t)+x(t)=0,x(t),t=-10..20,{[x(0)=2,D(x)(0)=-3]},stepsize=.2,linecolor=black);**

Note that the two initial conditions appear inside a
*single pair of brackets*
(because they determine
*only one solution curve*
). The stepsize specification was necessary, due to the oscillatory nature of the graph and the length of the plot interval; the default stepsize would have been 30/20 = 1.5 ! (Check out what you would get without stepsize specification.)

Basically, you have the same options as in the first-order case (except for those relating to direction fields, of course), and just as in the first-order case, several approximate solutions can be generated simultaneously, by specifying several sets of initial data in the form

{[x(t0)=x0,D(x)(t0)=v0],[x(t1)=x1,D(x)(t1)=v1],[x(t3)=x3,D(x)(t3)=v3],...} .

Note that, in analogy to the first-order case, the entire set of initial conditions is enclosed in curly braces, while the initial conditions for each individual solution curve appear in brackets. Example:

`> `
**DEplot(D(D(x))(t)+x(t)=0,x(t),t=-10..20,{[x(0)=2,D(x)(0)=-3],[x(0)=-3,D(x)(0)=2]},linecolor=black);**

Pretty lousy! Obviously, we should reduce the stepsize or use a more flexible solver. Let's do both and also add some perks:

`> `
**DEplot(D(D(x))(t)+x(t)=0,x(t),t=-10..20,{[x(0)=2,D(x)(0)=-3],[x(0)=-3,D(x)(0)=2]},
method=rkf45,stepsize=.2,linecolor=[blue,green],
title=`x(0)=2, x'(0)=-3 (blue); x(0)=-3, x'(0)=2 (green)`,
titlefont=[TIMES,ROMAN,12]);**