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. x ' = -k*x (the radioactive-decay equation), where the unknown function is x = x(t) and k is a parameter:

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

x(t) = _C1*exp(-k*t)

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

x(t) = c*exp(-k*t) .

The next command produces a general solution for the d.e. x '' + a^2*x = 0 (the elastic-spring equation with a^2 = k/m ), where the unknown function is again x = x(t) and a is a parameter:

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

x(t) = _C1*sin(a*t)+_C2*cos(a*t)

Note that the second derivative, x '', should be written as D(D(x)) 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. x ' = -k*x , x(0) = x0 :

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

x(t) = x0*exp(-k*t)

Next, we solve the second-order i.v.p. x '' + a^2*x = 0 , x(0) = x0 , x'(0) = v0 :

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

x(t) = v0*sin(a*t)/a+x0*cos(a*t)

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

> s:=rhs(%);

s := v0*sin(a*t)/a+x0*cos(a*t)

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, sol ) 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));

sol := x(t) = v0*sin(a*t)/a+x0*cos(a*t)

> s:=rhs(sol);

s := v0*sin(a*t)/a+x0*cos(a*t)

If we now specify the value of the parameter a and the initial values x0 and v0 , we can plot s just like any other function. Set, for example, a = 1 , x0 = 2 , and v0 = -3 :

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

a := 1

x0 := 2

v0 := -3

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

> s;

-3*sin(t)+2*cos(t)

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

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

[Maple Plot]

To obtain the solution and its graph for other values of the parameter a or other initial values x0 and v0 , 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 x ' = f(t,x) , in the rectangle t = a .. b , x = c .. d , 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 x ' = x^2 in the rectangle t = -3 .. 3 , x = -6 .. 6 :

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

[Maple Plot]

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

[Maple Plot]

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 x ' = x^2 , as before, along with the (approximate) graph of a particular solution, satisfying the initial condition x(0) = 1 :

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

[Maple Plot]

Now we will add a second solution curve, satisfying the initial condition x(0) = -1 :

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

[Maple Plot]

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 ( x -range) is optional; if omitted, DEplot will choose the x -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 x -range specification):

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

[Maple Plot]

Can you explain the poor outcome? (Hint: The exact solution with initial condition x(0) = 1 has a vertical asymptote at t = 1 .)

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);

[Maple Plot]

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`);

[Maple Plot]

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 ( b-a )/20, where a .. b is the horizontal range ( t -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 ( b-a )/20 . Let's try the following example, where we approximate the same solutions of x ' = x^2 as before, but with a larger vertical range ( -10 .. 10 ) 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);

[Maple Plot]

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, ( b-a )/20 with b = 3 and a = -3 ), 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);

[Maple Plot]

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.

dx/dt = 1-exp(t)*x^2 , x(0) = 0 ,

on the interval t = 0 .. 15 , 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);

[Maple Plot]

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);

[Maple Plot]

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 '' + x = 0 , x(0) = 2 , x'(0) = -3, whose symbolic solution is 2*cos(t)-3*sin(t) :

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

[Maple Plot]

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);

[Maple Plot]

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

[Maple Plot]