Lesson01.mw
DIFFERENTIAL EQUATIONS POWERTOOL
Lesson 1 -- Introduction to Differential Equations in Maple
Prof. Douglas B. Meade
Industrial Mathematics Institute
Department of Mathematics
University of South Carolina
Columbia, SC 29208
URL: http://www.math.sc.edu/~meade/
E-mail: meade@math.sc.edu
Copyright 2001 by Douglas B. Meade
All rights reserved
-------------------------------------------------------------------
Outline for Lesson 1
1.A Entering differential equations in Maple
1.A-1 First-order differential equations
1.A-2 Higher-order differential equations
1.A-3 Systems of differential equations
1.B Explicit solutions to differential equations
1.B-1 General solutions
1.B-2 Solutions to initial value problems (IVPs)
1.C Numerical Solution to an IVP
1.D Graphical Solution to an IVP
1.E Using
odeadvisor
to classify an ODE
Initialization
> |
interface( rtablesize=30 ); |
Warning, the name changecoords has been redefined
1.A Entering differential equations in Maple
A differential equation can be entered in Maple using any of the methods for constructing algebraic, transcendental, or any other
equation
in Maple. It is a good idea to assign each differential equation to a unique, and descriptive, Maple name. Such assignments are typically done using an
assignment statement
. For example, the ordinary differential equation (ODE)
would be entered into Maple as
> |
ode := diff( x(t), t ) = k * x(t); |
The
assign
command will not be of too much use in creating differential equations. Note that if
assign
is used to convert an equality into an assignment, the result of the assignment will not have an equal sign, i.e., it will not be an
equation
. When Maple expects an equation and receives an expression, the expression is assumed to be the left-hand side of an equation whose right-hand side is zero. For example, when working with Newton's second law of motion, one might write
> |
eq_of_motion := F = m*a; |
for the equation of motion, then set
, the acceleration, as
to obtain
> |
eq_of_motion2 := eval( eq_of_motion, a=diff(x(t),t,t) ); |
When the
assign
command is applied to eq_of_motion2,
the result is that the Maple name F has been assigned a value,
Note that F is only an expression, not an equation. In most instances, when Maple receives an expression when an equation is expected, it is assumed that the expression is the left-hand side of the equation and the right-hand side is zero. Thus, the following two commands are equivalent:
> |
dsolve( F, x(t) );
dsolve( F=0, x(t) ); |
1.A-1 First-order differential equations
A first-order differential equation is an equation that expresses a relationship between a function,
, its independent variable,
, and
, the first derivative. Other numeric or symbolic parameters can also appear in the equation.
The
diff
and
D
commands are used to represent derivatives in Maple, the first appearing as
> |
q1 := `x'` = diff (x(t), t ): |
and the second, as
The general form for a first-order differential equation is
.
In most cases it will be possible to solve for the derivative, i.e,
. While
dsolve
can generally work with ODEs specified in either form, the graphical command
DEplot
requires the explicit format (or its equivalent for a system).
1.A-2 Higher-order differential equations
A second-order derivative of a function
can be entered using
diff
as
> |
q2 := `x''` = diff( x(t), t,t ): |
The corresponding derivative expressed in terms of
D
is
If the right-hand side of this expression were entered explicitly, the syntax would be one of
or
The second form is functionally equivalent to the
composition
of D with itself, that is,
To represent higher-order derivatives, it is sometimes more convenient to use the dollar operator,
$
, instead of explicitly listing the independent variable in the
diff
command once for each derivative. Thus,
, the third derivative, would be entered as
> |
`x'''` := diff( x(t), t$3 ); |
The equivalent expressions in terms of the
D
command are
and
1.A-3 Systems of differential equations
The general form for a system of
first-order differential equations for the
functions
with independent variable
is
(
,
,
,
,
,
,
,
,
) = 0
(
,
,
,
,
,
,
,
,
) = 0
.
.
.
(
,
,
,
,
,
,
,
,
) = 0.
In most cases of interest, it will be possible to solve explicitly for each of the
first derivatives in terms of
and the
function values so that the ODEs can be put into the form
=
(
,
,
,
,
) for
= 1, 2,
,
.
While it is possible to express a system of differential equations in terms of Maple vectors, it is generally preferable to express the system as a list (or set) of
equations. For example, the second-order differential equation
> |
ode2 := diff( x(t), t,t ) + sin( x(t) ) = 0; |
can be represented as a system in terms of
and
, that is, as
which would appear in Maple as
> |
sys := [ diff( x(t), t ) = v(t), |
> |
diff( v(t), t ) = -sin(x(t)) ]; |
For a linear system, the set of explicit differential equations can be constructed from a coefficient matrix such as
The vector of unknowns would appear in Maple as
> |
X := < x[1](t), x[2](t) >; |
and the derivative of this vector of unknowns would appear as
> |
`X'` := map( diff, X, t ); |
The individual ODEs of the system can be generated in Maple via
> |
sys1 := eval( GenerateEquations( A, [_x1,_x2], `X'`), |
and then put into the more standard form
> |
sys2 := map( isolate, sys1, diff ); |
1.B Explicit solutions to differential equations
The
dsolve
command is used to obtain a solution to a differential equation. If initial and/or boundary conditions are specified, Maple attempts to find a particular solution to the specified initial or boundary value problem. Otherwise, the result is a general solution to the differential equation.
The terms general solution, and initial value problem (IVP) will be clarified below. The terms particular solution and boundary value problem (BVP) will be explained at appropriate places in the sequal.
1.B-1 General solutions
In its simplest form, the
dsolve
command accepts the differential equation and the unknown function. If the differential equation is
> |
ode := diff( x(t), t ) = k*x(t); |
and the unknown function is
then the solution is found in Maple with
> |
sol := dsolve( ode, fn ); |
The methods considered and executed by Maple to obtain a solution to an ODE or IVP can be observed by increasing the
infolevel
for
dsolve
to 3 :
> |
infolevel[dsolve] := 3: |
Then, re-executing the previous example yields :
`Methods for first order ODEs:`
`--- Trying classification methods ---`
`trying a quadrature`
`trying 1st order linear`
`<- 1st order linear successful`
To suppress the additional information, reset the
infolevel
for
dsolve
to 0 :
> |
infolevel[dsolve] := 0: |
Notes:
all occurrences of the dependent variable must appear with the independent variable, so the usage is
, not just
;
the general solution of this first-order ODE is a one parameter family of solutions where the default name of the parameter is _C1;
the result of
dsolve
is an equation; to access the right-hand side of this equation use the right-hand side command, namely,
rhs
:
applying the
assign
command to the result of
dsolve
does not convert the solution to a Maple function of the independent variable; to obtain the solution as a function of the independent variable, use
unapply
:
> |
X := unapply( rhs(sol), t ); |
Maple's unapply command converts an existing expression into a function. Since
is now a function, Maple can correctly interpret the symbols
as
The syntax for obtaining the general solution for a system of ODEs is similar. If the system is given by
> |
sys := diff( x(t), t ) = v(t), |
> |
diff( v(t), t ) = 2*x(t) - v(t); |
wherein the dependent variables are
the solution can be obtained via
> |
sol := dsolve( {sys}, {fns} ); |
Notes:
the first argument to dsolve must be either a single ODE or a set or list of ODEs; the second argument, the dependent function(s), follows a similar format
because the result returned by
dsolve
is a set of equations, use
subs
or
eval
to extract the individual components of the solution. (This approach can also be used in place of
rhs
for single equations. In fact, it is the more efficient strategy for extracting expressions from Maple computations.)
> |
X := subs( sol, x(t) );
V := eval( v(t), sol ); |
1.B-2 Solutions to initial value problems (IVPs)
An initial value problem is an ODE together with an appropriate set of initial conditions. To find the solution of an IVP with dsolve, the first argument must be a set consisting of the ODE and the initial conditions. If the ODE is
and the initial condition is
the dsolve solution is obtained via
> |
dsolve( {ode,ic}, fn ); |
An IVP for a system of ODEs requires one initial condition for each ODE. Thus, if the system is
and the initial conditions are
> |
sys_ic := x(0)=X0, v(0)=0; |
the dsolve solution is obtained via
> |
sys_soln := dsolve( {sys, sys_ic}, {fns} ); |
If the IVP contains a differential equation of order
, that is, the highest-ordered derivative in the ODE is an
th derivative, then there must be
initial conditions. These conditions will be prescribed values for the dependent variable and its first
derivatives all prescribed at a single point.
The following example is the second-order IVP equivalent to the first-order system just solved. The differential equation is
> |
ode2 := diff( x(t), t,t ) + diff( x(t), t ) - 2*x(t) = 0; |
and the two initial conditions are
> |
ic2 := x(0)=X0, D(x)(0)=0; |
The dsolve solution for
is then
> |
sol2 := dsolve( {ode2,ic2}, x(t) ); |
The
odetest
command can be used to verify the equivalence of the second-order ODE and the system of first-order ODEs. First, verify that the function
returned as part of the system's solution, namely,
> |
sys_solnX := op(select(has,sys_soln,x(t))); |
satisfies the second-order ODE
> |
odetest( sys_solnX, ode2 ); |
Alternatively, the function
satisfying the system of ODEs should satisfy
. Hence, we compute
, obtaining
> |
sol2v := v(t) = diff(rhs(sol2),t); |
The check that the functions
,
satisfy the system of ODEs is completed with the command
> |
odetest( {sol2,sol2v}, {sys} ); |
Note that Maple's
dsolve
command is not savvy enough to know what to do with a differential equation that does not explicitly identify the independent variable. For example, if the incorrect notation is used, dsolve returns an error, as we see with
> |
dsolve( D(x)=k*x, x(t) ); |
Error, (in ODEtools/info) Not an ODE w.r.t. x(t)
The correct syntax for this example would be
> |
dsolve( D(x)(t)=k*x(t), x(t) ); |
Maple supports different notation for an initial condition of the form
. The best notation to use in the context of differential equations is the
D
operator, the syntax for which is
The
D
operator can be replaced by the
operator by using the Maple command
> |
convert(D(x)(0)=2,diff); |
but it is not practical to use an expression of this type as an initial condition in the dsolve command. That it works is clear from
> |
dsolve( {diff(x(t),t)=k*x(t), convert(D(x)(0),diff)=1}, x(t) ); |
but most users prefer to use
diff
when creating a differential equation and
D
when creating initial (or boundary) conditions.
1.C Numerical solutions to an IVP
It is not possible to find an explicit solution for many IVPs. And, even when it is possible to find an explicit solution, the solution might be so complicated as to be of little or no practical use. In either event, Maple can be instructed to return a
numerical solution
to the IVP by including
type=numeric
as an optional argument to
dsolve
.
Note that values must be assigned to all parameters in the IVP before asking Maple for a numerical solution. Thus if the IVP is given by
> |
ode2;
ic3 := eval(ic2,X0=3); |
a numeric solution in Maple begins with
> |
nsol := dsolve( {ode2,ic3}, x(t), type=numeric ); |
Maple has created a procedure (subroutine) whose invocation, or call, will cause numbers to be generated. The procedure nsol is invoked at a value of the independent variable, here,
. For example, at
, the call would be
The valued returned by this procedure is a list of equations. The first equation specifies the value of the independent variable, subsequent equations display the value of each of the functions (including derivatives) for which an initial condition is specified. Note that each appearance of the dependent variable and its derivatives is
and not
. This allows values of
to be extracted with the substitution-device illustrated earlier, so that
could be referenced as
Again, note the use of the generic
to obtain the specific value of
.
The
odeplot
command, found in the
plots
package, can be used to produce one or more plots of the numeric solution to an IVP. In Figure 1.1 below, the graphs of the (approximate) solution (red) and its first derivative (green) are displayed on the interval [0, 1].
> |
odeplot( nsol, [[t,x(t)],[t,diff(x(t),t)]], 0..1, legend=[`x(t)`,`x'(t)`], title="Figure 1.1" ); |
A powerful feature of
odeplot
is its ability to plot functions of the dependent variable. For example, the graph of
on [0, 10], shown in Figure 1.2, provides strong evidence that
.
> |
odeplot( nsol, [t,x(t)-diff(x(t),t)], 0..10, title="Figure 1.2" ); |
The default numerical integration algorithm is a Fehlberg fourth-fifth Runge-Kutta method (
RKF45
). Both higher-order and classical methods are specified with the
method=
optional argument. For example, the approximate solution using Euler's method, with stepsize
, is obtained with
> |
Xeuler := dsolve( {ode2,ic3}, x(t), type=numeric, method=classical, stepsize=0.1 ); |
On the interval
, Figure 1.3, a plot comparing the RKF45 solution (in red) and Euler's method solution (in blue), is obtained with
> |
p1 := odeplot( nsol, [t,x(t)], 0..2, color=red ): |
> |
p2 := odeplot( Xeuler, [t,x(t)], 0..2, color=blue ): |
> |
display( [p1,p2], title="Figure 1.3" ); |
The dsolve command can also create a table of values for the solution of an IVP using any of the available methods. For example, the table of values of the solution (and its derivative) using the forward Euler method at the times in the list
> |
times := Array( [i/10$i=0..10, $2..10] ); |
is created with the command
> |
t1 := dsolve( {ode2,ic3}, x(t), type=numeric, method=classical[foreuler], output=Array([i/10$i=0..10, $2..10]) ); |
Note that this structure is a 2 x 1 matrix in which the (2,1) entry is itself a matrix. Thus, for example, the approximate value of the solution at
= 0.5 is
A more interesting use of a table of values is to compare the solutions for two different numerical methods (or for one method with different stepsize). To illustrate, the table computed using the Adams-Bashforth-Moulton predictor-corrector method is computed with
> |
t2 := dsolve( {ode2,ic3}, x(t), type=numeric, method=classical[abmoulton], output=Array([i/10$i=0..10, $2..10]) ): |
The columns of approximations from the forward Euler and predictor-corrector methods can be compared in a table created via
> |
header1 := <` t `| `Forward Euler`| `Predictor-Corrector`>: |
> |
header2 := <`---`| `-------------`| `-------------------`>: |
> |
< header1, header2, <Column(t1[2,1],1..2) | Column(t2[2,1],2) > >; |
1.D Graphical solutions to an IVP
One of the graphical tools for analyzing the solution of an IVP,
odeplot
, was discussed in Section 1.C . The
DEplot
command, from the
DEtools
package, provides graphical information directly from the differential equation(s). In it's simplest form,
DEplot
can create a plot of the solution to an IVP. If the IVP is
then DEplot graphs its solution in Figure 1.4 via the syntax
> |
DEplot( ode2, x(t), t=-1..2, [[ic3]], title="Figure 1.4" ); |
The extra set of square brackets around the initial condition are needed because
DEplot
facilitates the display of the solution to several different initial conditions in a single plot. For example, the plot of the solution for two additional sets of initial conditions can be added with one simple change to the
DEplot
command. If the additional initial conditions are entered as
> |
ic4 := x(0)=0, D(x)(0)=3; |
> |
ic5 := x(0)=-3, D(x)(0)=0; |
then graphs of solutions of three different IVPs posed for the same ODE appear in Figure 1.5 via the syatax
> |
DEplot( ode2, x(t), t=-1..2, [[ic3],[ic4],[ic5]], linecolor=[black,red,blue], title="Figure 1.5" ); |
Note how the individual solution curves are rendered in different colors by listing separate colors for the parameter linecolor.
To illustrate the capabilities of DEplot for a system of ODEs, consider the classical Lotka-Volterra model for a predator-prey interaction between a rabbit population,
, and a fox population,
(see also, Lesson 26). The ODEs for the system are
> |
sys_PP := diff(r(t),t)=r(t)*(2-f(t)), |
> |
diff(f(t),t)=.3*f(t)*(r(t)-1); |
and the dependent variables are
A direction field for this system is a graph of the
-plane in which are displayed arrows whose slopes are given by
and whose directions are determined by the increasing value of
, the independent variable. Thus, the arrows are tangent to solution curves and therefore give a sense how a collection of such curves would "flow."
The direction field for this system can be obtained with the command:
> |
DEplot( [sys_PP], [fns_PP], t=0..1, r=0..3, f=0..3, arrows=medium, |
> |
title="Direction Field for Predator-Prey System" ); |
Note that no initial conditions have been specified and that, because the system is autonomous (the independent variable does not appear explicitly in the ODEs), the range for the independent variable is immaterial to the appearance of the direction field (but is required by
DEplot
).
A phase portrait is a graph in the phase plane (here, the
-plane) that shows multiple solution solution curves, usually called orbits, trajectories, or paths. It is customary to include arrows showing, along these curves, the direction of increasing values of the independent variable.
Solutions that are to pass through the following initial points
> |
ic_PP := [r(0)=2.5,f(0)=i] $ i=1..5; |
are incorporated into a single graph with the command
> |
DEplot( [sys_PP], [fns_PP], t=0..30, [ic_PP], arrows=SMALL, |
> |
stepsize=0.1, linecolor=BLUE, |
> |
title=`Phase Portrait for Predator-Prey System` ); |
In the time-domain, solution curves for rabbit and fox populations passing through the first of these initial points can be obtained by specifying the optional argument scene. In Figure 1.6, the rabbit population
is shown in blue, whereas the fox population
is shown in red.
> |
pR := DEplot( [sys_PP], [fns_PP], t=0..30, [ic_PP[1]], |
> |
stepsize=0.1, linecolor=BLUE, scene=[t,r] ): |
> |
pF := DEplot( [sys_PP], [fns_PP], t=0..30, [ic_PP[1]], |
> |
stepsize=0.1, linecolor=RED, scene=[t,f] ): |
> |
display( [pR,pF], scaling=constrained, title="Figure 1.6" ); |
Additional optional arguments to
DEplot
allow for the specification of the numeical method to be used (the default is a fourth-order Runge-Kutta method,
method=classical[rk4]
), the color, number, and size of arrows in a direction field, and many other features beyond the needs for this course.
1.E Using
odeadvisor
to classify an ODE
The
odeadvisor
command, part of the
DEtools
package, provides a method for checking the classification of an ODE prior to obtaining a solution. The basic syntax for
odeadvisor
is very simple. If the differential equation is
> |
ode := diff( y(t), t ) = a*y(t); |
then its classification is obtained with
> |
odeadvisor( ode, y(t) ); |
This result is a little surprising. Most users, including the author himself, generally classify this ODE as either separable or first-order linear. When odeadvisor returns an unexpected, or unfamiliar, response it can be useful to include help as the final argument to odeadvisor, as shown by the following usage.
> |
odeadvisor( ode, y(t), help ); |
To check if an ODE is solvable by one or more specific methods, a list of methods can be included as an optional argument. Hence, to determine if the ODE is separable, use
> |
odeadvisor( ode, y(t), [separable] ); |
whereas to show that the ODE is not exact, use
> |
odeadvisor( ode, y(t), [exact] ); |
A (seemingly) simple modification to the ODE can significantly alter the nature, and hence, the classification of the ODE. For example, modifying the simple ODE above to
> |
ode2 := diff( y(t), t ) = a*y(t) + f(t); |
means its type changes to
a form that is no longer separable, as seen with
> |
odeadvisor( ode2, [separable] ); |
Second-, and higher-, order ODEs can be handled in a similar manner. For example, the second-order ODE
> |
ode3 := x^2*diff(y(x),x$2) + 2*x*diff(y(x),x) + y(x) = f(x); |
is classified as
> |
odeadvisor( ode3, y(x) ); |
The following usage will yield a help-page full of additional information.
> |
odeadvisor( ode3, help ); |
Unfortunately, odeadvisor is not able to provide any assistance for a system of ODEs. For example, if the system
> |
sys1 := { diff(v(t),t) + x(t) = 0, diff(x(t),t)=v(t) }; |
is given to odeadvisor, an error results, as seen from
> |
odeadvisor( sys1, {x(t),v(t)} ); |
Error, (in odeadvisor) invalid input: `ODEtools/odeadv` expects its 1st argument, ODE, to be of type `ODEtools/ODE`, but received {diff(x(t), t) = v(t), (diff(v(t), t))+x(t) = 0}
Of course, if the system is converted to an equivalent higher-order ODE, odeadvisor might be able to provide some useful information. When the system is written as the second-order ODE
> |
ode4 := diff(x(t),t$2) + x(t) = 0; |
Maple returns
It is not a surprise to have odeadvisor return results that are beyond the scope of a first course, or any undergraduate course, in differential equations. For example, the ODE
> |
ode5 := (2*y(x)-x)*diff(y(x),x)-y(x)-2*x; |
provokes the response
This means that we will expect Maple to be able to solve many ODEs and IVPs that we won't be able to solve by hand at the end of this course. It also means that when we can solve an equation by hand, our results may appear different than Maple's because of the different methods used to obtain the solution. Recall that
infolevel[dsolve]:=3:
causes Maple to display information about the methods used by the
dsolve
command;
infolevel[dsolve]:=0:
turns off these messages.
[Back to ODE Powertool Table of Contents]