DIFFERENTIAL EQUATIONS POWERTOOL

Unit 1 -- Introduction to Differential Equations in Maple

Outline for Unit 1

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

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 .

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

> eq_of_motion := F = m*a;

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

>

A first-order differential equation is an equation that expresses a relationship between a function, x, its independent variable, , and the derivative, . (Note that Maple expresses all derivatives as partial derivatives. This is a fact of life that we must accept.) Other numeric or symbolic parameters can appear in the equation. The diff and D commands are used to represent derivatives in Maple.

> q1 := `x'` = diff (x(t), t ):

> q1;

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

>

A second-order derivative of a function x 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.

> `x'''` := diff( x(t), t$3 );

>

The equivalent expressions in terms of the D command are

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

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

>

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 explicitly solve for each of the first derivatives in terms of and the function values. That is,

= ( , , , ..., ) 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

> 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 the coefficient matrix.

> A := matrix( 2, 2, [[a,b],[c,d]] );

> X := vector( 2, [x[1](t), x[2](t)] );

> `X'` := map( diff, X, t );

> sys1 := geneqns( A, X, `X'`);

> sys2 := map( isolate, sys1, diff );

>

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 simplest form of the dsolve command requires the differential equation and the unknown function:

> ode := diff( x(t), t ) = k*x(t);

> fn := x(t);

> 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 to isolate the derivative dx/dt...

Successful isolation of dx/dt

-> 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 unknown function in the ODE include the independent variable, , not just the function name, ;
• since this is a general solution to a first-order equation, the general solution is a one parameter family of solutions, the default name of the parameter is _C1 ;
• the result of dsolve is an equation, to access the right-hand side of this definition use 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 );

> X(0);

>

The syntax for obtaining the general solution for a system of ODEs is similar:

> sys := diff( x(t), t ) = v(t),

> diff( v(t), t ) = 2*x(t) - v(t);

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

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

> X := subs( sol, x(t) );

> V := eval( v(t), sol );

>

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:

> ode;

> ic := x(0)=A;

> dsolve( {ode,ic}, fn );

>

An IVP for a system of ODEs requires one initial condition for each ODE.

> sys;

> sys_ic := x(0)=X0, v(0)=0;

> sys_soln := dsolve( {sys, sys_ic}, {fns} );

>

For an IVP for an ODE of order n requires n initial conditions, i.e., values for the unknown function and its first (n-1) derivatives at a single point. The following example is the second-order IVP is equivalent to the the first-order system just considered.

> ode2 := diff( x(t), t,t ) + diff( x(t), t ) - 2*x(t) = 0;

> ic2 := x(0)=X0, D(x)(0)=0;

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

> 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 is computed to be

> sol2v := v(t) = diff(rhs(sol2),t);

>

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

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

Error, (in ODEtools/info) Not an ODE w.r.t. x(t)

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

>

> convert(D(x)(0)=2,diff);

This is just one of many reasons why most users prefer to use diff when creating a differential equation.

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

Error, (in ODEtools/info) Not an ODE w.r.t. x(t)

>

The D operator is generally used to specify derivative initial and/or boundary conditions. For example,

> dsolve( {diff(x(t),t)=k*x(t), D(x)(0)=1}, x(t) );

>

While the derivative condition can be converted to an equivalent expression in terms of diff,

> convert(D(x)(0)=1,diff);

>

it is not practical to use this expression as an initial condition in dsolve :

> dsolve( {diff(x(t),t)=k*x(t), convert(D(x)(0),diff)=1}, x(t) );

>

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 to be of little or no practical use. When such situations are encountered, 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.

> ode2;

> ic3 := eval(ic2,X0=3);

> nsol := dsolve( {ode2,ic3}, x(t), type=numeric );

>

The procedure requires a value of the independent variable as its argument. 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.

> nsol(1);

>

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 the first example, 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)`] );

>

A powerful feature of odeplot is the ability to plot expressions involving the solution of the IVP. Here, the plot of on [0,10] provides strong evidence that :

> odeplot( nsol, [t,x(t)-diff(x(t),t)], 0..10 );

>

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

>

A plot comparing the RKF45 and Euler's method solutions on the interval [0,2] 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] );

>

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

> 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], value=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 t=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], value=array([i/10$i=0..10, $2..10]) ):

>

The column corresponding to the solution, x(t), from the predictor-corrector method can be compared to the corresponding column for the forward Euler method using

> header1 := vector([` t `, `Forward Euler`, `Predictor-Corrector`]):

> header2 := vector([`---`, `-------------`, `-------------------`]):

> stackmatrix( header1, header2, augment( col(t1[2,1],1..2), col(t2[2,1],2) ));

>

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:

> ode2;ic3;

> DEplot( ode2, x(t), t=-1..2, [[ic3]] );

>

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:

> ic4 := x(0)=0, D(x)(0)=3;

> ic5 := x(0)=-3, D(x)(0)=0;

> DEplot( ode2, x(t), t=-1..2, [[ic3],[ic4],[ic5]] );

>

The odeadvisor command, part of the DEtools package, provides a method to check the classification of an ODE prior to obtaining a solution. The basic syntax for odeadvisor is very simple.

> ode := diff( y(t), t ) = a*y(t);

> odeadvisor( ode, y(t) );

>

This result is a little surprising. Most users, including myself, 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 .

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

> odeadvisor( ode, y(t), [separable] );

> odeadvisor( ode, y(t), [exact] );

>

A (seemingly) simple modification to the ODE can significantly alter the classification of the ODE. For example,

> ode2 := diff( y(t), t ) = a*y(t) + f(t):

> odeadvisor( ode2 );

> odeadvisor( ode2, [separable] );

>

Second-, and higher-, order ODEs can be handled in a similar manner.

> ode3 := x^2*diff(y(x),x$2) + 2*x*diff(y(x),x) + y(x) = f(x);

> odeadvisor( ode3, help );

>

Unfortunately, odeadvisor is not able to provide any assistance for a system of ODEs.

> sys1 := { diff(v(t),t) + x(t) = 0, diff(x(t),t)=v(t) };

> odeadvisor( sys1, {x(t),v(t)} );

Error, (in odeadvisor) ODEtools/odeadv expects its 1st argument, ODE, to be of type ODEtools/ODE, but received {diff(v(t),t)+x(t) = 0, diff(x(t),t) = v(t)}

>

Of course, if the system is converted to an equivalent higher-order ODE, odeadvisor might be able to provide some useful information.

> ode4 := diff(x(t),t$2) + x(t) = 0;

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

> ode5 := (2*y(x)-x)*diff(y(x),x)-y(x)-2*x;

> 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 from 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.

>