DIFFERENTIAL EQUATIONS POWERTOOLLesson 1 -- Introduction to Differential Equations in MapleProf. Douglas B. MeadeIndustrial Mathematics InstituteDepartment of MathematicsUniversity of South CarolinaColumbia, SC 29208 URL: http://www.math.sc.edu/~meade/E-mail: meade@math.sc.edu Copyright \251 2001 by Douglas B. MeadeAll rights reserved-------------------------------------------------------------------

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 equations1.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 IVP1.D Graphical Solution to an IVP1.E Using odeadvisor to classify an ODE

restart;with( DEtools ):with( plots ):with( LinearAlgebra ):interface( rtablesize=30 );

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)NiMvKiYlI2R4RyIiIiUjZHRHISIiKiYlImtHRiYlInhHRiY= would be entered into Maple asode := 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 writeeq_of_motion := F = m*a;for the equation of motion, then set NiMlImFH, the acceleration, as NiMtJSR4JydHNiMlInRH to obtaineq_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, NiMlInhH, its independent variable, NiMlInRH, and NiMtJSVkaWZmRzYkLSUieEc2IyUidEdGKQ==, 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 asq1 := `x'` = diff (x(t), t ):q1;and the second, asconvert( q1, D );The general form for a first-order differential equation is NiMvLSUiRkc2JSUidEctJSJ4RzYjRictJSVkaWZmRzYkRihGJyIiIQ==. In most cases it will be possible to solve for the derivative, i.e, NiMvLSUlZGlmZkc2JC0lInhHNiMlInRHRiotJSJHRzYkRipGJw==. 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 NiMlInhH can be entered using diff asq2 := `x''` = diff( x(t), t,t ):q2;The corresponding derivative expressed in terms of D isconvert( 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, NiMlJXgnJydH, 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);

The general form for a system of NiMlIm1H first-order differential equations for the NiMlIm5H functions NiYmJSJ4RzYjIiIiJkYkNiMiIiMlJC4uLkcmRiQ2IyUibkc= with independent variable NiMlInRH isNiMmJSJGRzYjIiIi(NiMlInRH, NiMtJiUieEc2IyIiIjYjJSJ0Rw==, NiMtJiUieEc2IyIiIzYjJSJ0Rw==, NiMlJC4uLkc=, NiMtJiUieEc2IyUibkc2IyUidEc=, NiMtJSVkaWZmRzYkLSYlInhHNiMiIiI2IyUidEdGLA==, NiMtJSVkaWZmRzYkLSYlInhHNiMiIiM2IyUidEdGLA==, NiMlJC4uLkc=, NiMtJSVkaWZmRzYkLSYlInhHNiMlIm5HNiMlInRHRiw= ) = 0NiMmJSJGRzYjIiIj(NiMlInRH, NiMtJiUieEc2IyIiIjYjJSJ0Rw==, NiMtJiUieEc2IyIiIzYjJSJ0Rw==, NiMlJC4uLkc=, NiMtJiUieEc2IyUibkc2IyUidEc=, NiMtJSVkaWZmRzYkLSYlInhHNiMiIiI2IyUidEdGLA==, NiMtJSVkaWZmRzYkLSYlInhHNiMiIiM2IyUidEdGLA==, NiMlJC4uLkc=, NiMtJSVkaWZmRzYkLSYlInhHNiMlIm5HNiMlInRHRiw= ) = 0...NiMmJSJGRzYjJSJtRw==(NiMlInRH, NiMtJiUieEc2IyIiIjYjJSJ0Rw==, NiMtJiUieEc2IyIiIzYjJSJ0Rw==, NiMlJC4uLkc=, NiMtJiUieEc2IyUibkc2IyUidEc=, NiMtJSVkaWZmRzYkLSYlInhHNiMiIiI2IyUidEdGLA==, NiMtJSVkaWZmRzYkLSYlInhHNiMiIiM2IyUidEdGLA==, NiMlJC4uLkc=, NiMtJSVkaWZmRzYkLSYlInhHNiMlIm5HNiMlInRHRiw= ) = 0.In most cases of interest, it will be possible to solve explicitly for each of the NiMlIm5H first derivatives in terms of NiMlInRH and the NiMlIm5H function values so that the ODEs can be put into the form NiMtJSVkaWZmRzYkLSYlInhHNiMlImlHNiMlInRHRiw= = NiMmJSJHRzYjJSJpRw==(NiMlInRH, NiMtJiUieEc2IyIiIjYjJSJ0Rw==, NiMtJiUieEc2IyIiIzYjJSJ0Rw==, NiMlJC4uLkc=, NiMtJiUieEc2IyUibkc2IyUidEc= ) for NiMlImlH = 1, 2, NiMlJC4uLkc=, NiMlIm5H.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 NiMlIm1H equations. For example, the second-order differential equationode2 := diff( x(t), t,t ) + sin( x(t) ) = 0;can be represented as a system in terms of NiMtJSJ4RzYjJSJ0Rw== and NiMvLSUidkc2IyUidEctJSVkaWZmRzYkLSUieEdGJkYn, that is, asNiMvLSUjeCdHNiMlInRHLSUidkdGJg== NiMvLSUjdidHNiMlInRHLCQtJSRzaW5HNiMtJSJ4R0YmISIi which would appear in Maple assys := [ 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 asA := < <a|b>, <c|d> >;The vector of unknowns would appear in Maple asX := < 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 viasys1 := eval( GenerateEquations( A, [_x1,_x2], `X'`), {_x1=X[1],_x2=X[2]} );and then put into the more standard formsys2 := 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 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.

In its simplest form, the dsolve command accepts the differential equation and the unknown function. If the differential equation isode := diff( x(t), t ) = k*x(t);and the unknown function isfn := x(t);then the solution is found in Maple withsol := 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 );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 NiMtJSJ4RzYjJSJ0Rw==, not just NiMlInhH;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 NiMlIlhH is now a function, Maple can correctly interpret the symbols NiMtJSJYRzYjIiIh asX(0);The syntax for obtaining the general solution for a system of ODEs is similar. If the system is given bysys := diff( x(t), t ) = v(t), diff( v(t), t ) = 2*x(t) - v(t);wherein the dependent variables arefns := x(t), v(t);the solution can be obtained viasol := 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 formatbecause 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 );

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 isunassign('A'); ode;and the initial condition isic := x(0)=A;the dsolve solution is obtained viadsolve( {ode,ic}, fn );An IVP for a system of ODEs requires one initial condition for each ODE. Thus, if the system issys;and the initial conditions aresys_ic := x(0)=X0, v(0)=0;the dsolve solution is obtained viasys_soln := dsolve( {sys, sys_ic}, {fns} );If the IVP contains a differential equation of order NiMlIm5H, that is, the highest-ordered derivative in the ODE is an NiMlIm5Hth derivative, then there must be NiMlIm5H initial conditions. These conditions will be prescribed values for the dependent variable and its first NiMsJiUibkciIiJGJSEiIg== 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 isode2 := diff( x(t), t,t ) + diff( x(t), t ) - 2*x(t) = 0;and the two initial conditions areic2 := x(0)=X0, D(x)(0)=0;The dsolve solution for NiMtJSJ4RzYjJSJ0Rw== is thensol2 := 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 NiMtJSJ4RzYjJSJ0Rw== returned as part of the system's solution, namely,sys_solnX := op(select(has,sys_soln,x(t)));satisfies the second-order ODEodetest( sys_solnX, ode2 );Alternatively, the function NiMtJSJ2RzYjJSJ0Rw== satisfying the system of ODEs should satisfy NiMvLSUidkc2IyUidEctJSN4J0dGJg==. Hence, we compute NiMtJSN4J0c2IyUidEc=, obtainingsol2v := v(t) = diff(rhs(sol2),t);The check that the functions NiMtJSJ4RzYjJSJ0Rw==, NiMtJSJ2RzYjJSJ0Rw== satisfy the system of ODEs is completed with the commandodetest( {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 withdsolve( D(x)=k*x, x(t) );The correct syntax for this example would bedsolve( D(x)(t)=k*x(t), x(t) );Maple supports different notation for an initial condition of the form NiMvLSUjeCdHNiMiIiEiIiM=. The best notation to use in the context of differential equations is the D operator, the syntax for which isD(x)(0) = 2;The D operator can be replaced by the NiMqJiUiZEciIiIlI2R0RyEiIg== operator by using the Maple commandconvert(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.

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 byode2; ic3 := eval(ic2,X0=3);a numeric solution in Maple begins withnsol := 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, NiMlInRH. For example, at NiMvJSJ0RyIiIg==, the call would bensol(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 NiMtJSJ4RzYjJSJ0Rw== and not NiMtJSJ4RzYjIiIi. This allows values of NiMlInhH to be extracted with the substitution-device illustrated earlier, so that NiMtJSJ4RzYjIiIi could be referenced aseval(x(t), nsol(1));Again, note the use of the generic NiMtJSJ4RzYjJSJ0Rw== to obtain the specific value of NiMtJSJ4RzYjIiIi.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 NiMsJi0lInhHNiMlInRHIiIiLSUlZGlmZkc2JEYkRichIiI= on [0, 10], shown in Figure 1.2, provides strong evidence that NiMvLSUmbGltaXRHNiQsJi0lInhHNiMlInRHIiIiLSUlZGlmZkc2JEYoRishIiIvRislKWluZmluaXR5RyIiIQ==.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 NiMvJSJoRy0lJkZsb2F0RzYkIiIiISIi, is obtained withXeuler := dsolve( {ode2,ic3}, x(t), type=numeric, method=classical, stepsize=0.1 );On the interval NiM3JCIiISIiIw==, Figure 1.3, a plot comparing the RKF45 solution (in red) and Euler's method solution (in blue), is obtained withp1 := 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 listtimes := Array( [i/10$i=0..10, $2..10] );is created with the commandt1 := 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 NiMlInRH = 0.5 ist1[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 witht2 := 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 viaheader1 := <` t `| `Forward Euler`| `Predictor-Corrector`>:header2 := <`---`| `-------------`| `-------------------`>:< header1, header2, <Column(t1[2,1],1..2) | Column(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. If the IVP isode2;ic3;then DEplot graphs its solution in Figure 1.4 via the syntaxDEplot( 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 asic4 := 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 syataxDEplot( 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, NiMtJSJyRzYjJSJ0Rw==, and a fox population, NiMtJSJmRzYjJSJ0Rw== (see also, Lesson 26). The ODEs for the system aresys_PP := diff(r(t),t)=r(t)*(2-f(t)), diff(f(t),t)=.3*f(t)*(r(t)-1);and the dependent variables arefns_PP := r(t), f(t);A direction field for this system is a graph of the NiMlI3JmRw==-plane in which are displayed arrows whose slopes are given by NiMvKiYlI2RmRyIiIiUjZHJHISIiKiwtJSZGbG9hdEc2JCIiJEYoRiYlImZHRiYsJiUickdGJkYmRihGJkYwRigsJiIiI0YmRi5GKEYo and whose directions are determined by the increasing value of NiMlInRH, 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 NiMlI3JmRw==-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 pointsic_PP := [r(0)=2.5,f(0)=i] $ i=1..5;are incorporated into a single graph with the commandDEplot( [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 NiMtJSJyRzYjJSJ0Rw== is shown in blue, whereas the fox population NiMtJSJmRzYjJSJ0Rw== 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.

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 isode := diff( y(t), t ) = a*y(t);then its classification is obtained withodeadvisor( 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, useodeadvisor( ode, y(t), [separable] );whereas to show that the ODE is not exact, useodeadvisor( 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 toode2 := 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 withodeadvisor( ode2, [separable] );Second-, and higher-, order ODEs can be handled in a similar manner. For example, the second-order ODEode3 := x^2*diff(y(x),x$2) + 2*x*diff(y(x),x) + y(x) = f(x);is classified asodeadvisor( 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)} );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 ODEode4 := diff(x(t),t$2) + x(t) = 0;Maple returnsodeadvisor( 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 ODEode5 := (2*y(x)-x)*diff(y(x),x)-y(x)-2*x;provokes the responseodeadvisor( 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.

[Back to ODE Powertool Table of Contents]