Application Center - Maplesoft

Lesson 1: Introduction to Differential Equations in Maple

You can switch back to the summary page by clicking here.

Lesson01.mw

DIFFERENTIAL EQUATIONS POWERTOOL

Lesson 1 -- Introduction to Differential Equations in Maple

Columbia, SC 29208

-------------------------------------------------------------------

 >

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

 > restart;

 > with( DEtools ):

 > with( plots ):

 > with( LinearAlgebra ):

 > 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,

 > assign(eq_of_motion2);

 >

the result is that the Maple name F has been assigned a value,

 > F;

 >

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

 > q1;

 >

and the second, as

 > convert( q1, D );

 >

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

 > q2;

 >

The corresponding derivative expressed in terms of D is

 > convert( q2, D );

 >

If the right-hand side of this expression were entered explicitly, the syntax would be one of

 > (D@D)(x)(t);

 >

or

 > (D@@2)(x)(t);

 >

The second form is functionally equivalent to the composition of D with itself, that is,

 > D( D(x) )(t);

 >

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

 > (D@D@D)(x)(t);

 >

and

 > (D@@3)(x)(t);

 >

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

 > A := < , >;

 >

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'`),

 > {_x1=X[1],_x2=X[2]} );

 >

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

 > fn  := x(t);

 >

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 :

 > dsolve( ode, fn );

`Methods for first order ODEs:`
`--- Trying classification methods ---`

`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 :

 > X := rhs(sol);

 >

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

 > X(0);

 >

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

 > fns := x(t), v(t);

 >

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

 > unassign('A'); ode;

 >

and the initial condition is

 > ic := x(0)=A;

 >

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

 > sys;

 >

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

 > D(x)(0) = 2;

 >

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

 > nsol(1);

 >

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

 > eval(x(t), nsol(1));

 >

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

 > t1[2,1][6,2];

 >

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 := <`---`| `-------------`| `-------------------`>:

 >

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

 > ode2;ic3;

 >

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

 > fns_PP := r(t), f(t);

 >

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

 > odeadvisor( ode2 );

 >

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

 > odeadvisor( ode4 );

 >

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

 > odeadvisor( ode5 );

 >

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.

 >