odes.html

> restart;

>

Differential Equations

In this worksheet we demonstrate basic techniques of working with ordinary differential equations , (ODE's), using Maple. An ODE is a differential equation in which the unknown function is a function of one variable. We shall stick with first order ODE's which contain only the first derivative of the unknown function. Moreover, we shall consider only first order ODE's in the normal form , that is in the form:

where the left hand side is simply the derivative of the unknown function, in this case, and the right hand side is a given function of x and y. You know that for such ODE's, under reasonable assumptions about , there exists a unique solution to any initial value problem, ( IVP ), of the form

, ,

where , are given.

>

Solving and Plotting ODE's

Basic tools for solving and plotting ODE's are contained in the packages " plots " and " DEtools ". We begin with loading these packages. Please, don't forget to click on the two commands below.

> with(plots):

> with(DEtools):

>

We are familiar with the package " plots ". If you are curious about the content of " DEtools ", replace the colon at the end of the command with a semicolon and click on it again.

Example 1. (a) Find the general solution to the ODE: .

(b) Solve the following two initial value problems:

, ,

and

, .

>

In order to simplify many commands below let's first label our ODE:

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

>

Observe that correct syntax. The derivative is entered using the " diff " command. The command " D(y)(x) " could be used as well. Note that " y " has to be entered as " y(x) ". The main command for solving ODE's is " dsolve ".

> dsolve(ODE1,y(x));

>

Maple returned the general solution. " _C1 " denotes, of course, an arbitrary constant. Instead of using the name "ODE1" you could have entered the differential equation "diff(y(x),x)=-2*x*y(x)" directly into the " dsolve " command. Maple can handle initial value problems, as well. The proper syntax looks as follows.

> dsolve({ODE1,y(0)=2},y(x)); dsolve({ODE1,y(0)=1/2},y(x));

>

Maple can plot slope fields, as well as slope fields together with particular solutions. Proper commands are " dfieldplot " and "DEplot ", both contained in the " DEtools " package. Let's see how they work. Pay attention to the syntax.

> dfieldplot(ODE1,[y(x)],x=-2..2,y=-2..2,color=blue,scaling=constrained,arrows=LINE,dirgrid=[30,30]);

>

Maple plotted the slope field for our equation. All the options under the " dfieldplot " command regarding color, appearance of the arrows, scaling and dirgrid are, of course, optional. You can play with them and see what will happen. " dirgrid " tells Maple how dense you want the field of slopes to be. The default setting is dirgrid[20,20] and it tends to be a little rough. On the other hand, a finer grid may take more time to compute.

Remark. Whenever plotting field of slopes, you should use the "scaling = constrained" option. Otherwise, the pictures may appear misleading, as slopes become distorted.

The command " DEplot " plots the slope field together with particular solutions.

> DEplot(ODE1,y(x),-2..2,[[y(0)=2],[y(0)=1/2]], linecolor=magenta,color=blue,scaling=constrained,arrows=LINE);

>

As we expected, the two particular solutions are bell-shaped curves. If you do not want the slope field plotted with particular solutions, you add an option " arrows=NONE " under the " DEplot " command.

>

Example 2. (a) Try to find the general solution to the ODE: .

(b) Solve the IVP:

, .

Let's see how Maple handles this ODE.

> ODE2:=diff(y(x),x)=cos(x*y(x));

> dsolve(ODE2,y(x));

>

Maple returned no output! That means Maple is unable to solve the equation. If you are curious what steps Maple went through to find a solution before failing to do so, you can ask to see the steps using the command " infolevel ". The levels of information that you can request range from 0 to 5. The higher number, the more information Maple will reveal.

> infolevel[dsolve]:=3:

> dsolve(ODE2,y(x));

Methods for first order ODEs:

trying to isolate the derivative dy/dx

successful isolation of dy/dx

 -> trying classification methods

trying a quadrature

trying 1st order linear

trying Bernoulli

trying separable

trying inverse linear

trying simple symmetries for implicit equations

trying homogeneous types:

trying Chini

trying exact

 -> trying 2nd set of classification methods

trying Riccati

trying Abel

 -> trying Lie symmetry methods

>

As you see, Maple tried to match the equation with the types of first order ODE's that it knows how to solve in a closed form, and failed. Out of curiosity, let's see how Maple solved the previous equation, ODE1 , with which it was successful.

> dsolve(ODE1,y(x));

Methods for first order ODEs:

trying to isolate the derivative dy/dx

successful isolation of dy/dx

 -> trying classification methods

trying a quadrature

trying 1st order linear

1st order linear successful

>

As you see, Maple solved the equation by the first successful method, that is, as a linear equation in y(x). You don't know that method yet. You could, however, solve the equation by hand as a separable equation. Besides the common type equations listed above, Maple is familiar with more sophisticated techniques of solving ODE's involving power series expansion, Laplace transform and others. You have to tell Maple if you want it to apply those techniques. We are not going to do so at this point. Let's get back to the normal infolevel for " dsolve " command.

> infolevel[dsolve]:=0:

>

Maple couldn't find the general solution of the equation and neither can we. However, the solution to the IVP of (b), can be found using numerical methods . We describe them in the next section.

Solving Initial Value Problems Numerically Using Maple

Maple is programmed for the whole array of sophisticated numerical methods for solving ODE's. Let's find the numeric solution to the IVP of Example 2. (b). The proper syntax for evoking Maple's numerical solver is " dsolve(....,numeric) ". You should always label the resulting solution. For example as " soln ".

> soln:=dsolve({ODE2,y(0)=3},y(x),numeric);

>

The output looks like a failure. It is not. Maple reports its numerical solution as an algorithm that allows us to calculate values of the solution at any point we want, as well as plot it. The expression " rkf45 " stands for the Fehlberg fourth-fifth order Runge-Kutta method, which is a well-known numerical scheme for solving ODE's. Maple uses it as a default option. You can guess that the rkf45 algorithm is much more advanced that the Euler method that we have learned.

We can find values of our numerical solution " soln " at any single point we want, or at a list of points using the following syntax:

> [soln(1),soln(1.5),soln(2),soln(2.5),soln(3),soln(3.5),soln(4),soln(4.5),soln(5)];

>

Maple displays the consecutive values of x together with the corresponding values of the solution y(x). We can also plot a numerical solution to an IVP using the " odeplot " command. This command is contained in the " plots " package, which we have already loaded. The command can be used only with numerical solutions.

> odeplot(soln,[x,y(x)],0..5,color=magenta,thickness=2,scaling=constrained);

>

It may be interesting to see how this numerical solution relates to the slope field of the equation . To see both together you could use the " DEplot " command. You can also use the " display " command from the " plots " package.

> p1:=odeplot(soln,[x,y(x)],0..5,color=magenta,thickness=2,scaling=constrained):

> p2:=dfieldplot(ODE2,[y(x)],x=0..5,y=0..5,color=blue,arrows=LINE,scaling=constrained,dirgrid=[30,30]):

> display([p1,p2]);

>

An Applied Example: Mr. Jones' Date

Example 3. Mr. Jones has invited his date for dinner for 6 p.m. He plans to serve baked chicken. At 5:30 p.m. he realizes that he hasn't started baking the chicken. In a panic, Mr. Jones forgets to preheat the oven. He puts cold chicken from the fridge, at F, into a cold oven and turns the oven on. The temperature of the oven is rising according to the function

where t is the time, in minutes, after the oven is turned on. The internal temperature of the chicken, , in degrees Fahrenheit, t minutes after the oven is turned on, changes according to Newton's Law of Cooling:

The chicken is done when its internal temperature reaches F. At what time will Mr. Jones be able to serve the dinner?

>

To find the function , we have to solve the following IVP:

, .

Then we can determine when . Observe that this equation is different from the other ones we have seen in connection with Newton's Law of Cooling because in our present example we are not assuming that the temperature of the oven is constant. As a result, the differential equation is no longer separable, we will have trouble solving it by hand. Let's hope Maple can help us.

> ODE3:=diff(H(t),t)=-(7/1000)*(H(t)-350+280*exp(-(7/10)*t));

Maple somewhat simplified the equation. Let's solve the IVP.

> dsolve({ODE3,H(0)=40},H(t));

>

We have replaced all decimals by fractions in our ODE, as in earlier releases of Maple the " dsolve " command does not work if an equation contains decimals.

IMPORTANT NOTE: It is very important to realize that the output of the "dsolve" command may seem like a function, but, as far as Maple is concerned, H(t) is neither a function nor an expression. H(t) has not been properly defined as either one. In order to further process a solution to a given ODE, you have to make it into an expression or a function.

To obtain your solution as an expression, say exH , you can use the following syntax:

> dsolve({ODE3,H(0)=40},H(t)); exH:=rhs(%);

>

" rhs " stands for the " right hand side ". Recall, that " % " stands for the last output. Now, exH is an expression in terms of t and you can use commands, like " plot " and others with it. If you prefer to have your solution as a function, and you don't want to type long formulas, you have to use the command with a strange name " unapply " . The command turns an expression into a function as follows. We shall turn the expression exH in terms of t into a function fH of t.

> fH:=unapply(exH,t);

>

The logic behind the name "unapply" is that an expression in terms of t may be thought of as a result of applying a certain function to a given t. Hence, to turn an expression back into a function we have to "unapply".

Let's plot the function fH to see what is happening to Mr. Jones' chicken.

> plot(fH(t),t=0..90);

>

We see the the temperature of the chicken will reach 175 degrees somewhere between t = 60 and t = 90.

> fsolve(fH(t)=175,t,60..90);

>

Well, the chicken will be ready about 83 minutes after 5:30 p.m., that is, about 6:50 p.m. Mr. Jones has to hope that his date doesn't arrive too hungry!

>

Euler's Method

Maple can be used to illustrate Euler's method that we have learned. One can say that Euler's method is to numerical methods for ODE's what Left Sums are to numerical integration. It is a rudimetary and not very efficient numerical method, which, however, is easy to understand and very illuminating. Maple has Euler's method built in as one of the options, but in order to see all the steps, we shall program it from scratch.

Example 4. Construct approximate solutions for x from 0 to 5 to the initial value problem

using Euler's method with the three different step sizes = 0.5, 0.25, 0.125. Compare your solutions with Maple's solution to the IVP.

>

To execute Euler's method, we could go three times step-by-step through ten steps. We could tell Maple in one command to do the ten steps for our particular IVP using the so-called loop structure. Instead, we shall program a simple procedure that calculates Euler approximations for any given ODE with a right hand side f(x,y), any initial conditions , and any number od steps n. The loop structure " do...od " is a part of the procedure. This is a simple example of Maple programming. We will program a procedure called "Eulermethod" which, given f(x,y), initial conditions (x0,y0), step size h, and number of steps n, will calculate the steps of Euler's method and display the resulting list of pairs [ [x[0],y[0]], [x[1],y[1]], ...,[x[n],y[n] ]. To clarify the structure of our little program we will comment on each step.

> Eulermethod:=proc(f,x0,y0,h,n) local i,x,y;
x[0]:=x0; y[0]:=y0; first pair of the list is given by the initial condition
for i from 1 to n we tell Maple to perform n steps
do mark beginning of the loop commands
y[i]:=y[i-1]+evalf(f(x[i-1],y[i-1]))*h; get the next y value
x[i]:=x[i-1]+h; get the next x value
od; mark the end of the loop commands
[seq([x[i-1],y[i-1]],i=1..n+1)]; form the final list of pairs
end:

>

Let's apply the procedure to our function .

> f:=(x,y)->cos(x*y);

>

Lat's use Eulermethod with step sizes h = 0.5 =1/2 , h = 0.25=1/4, h=0.125=1/8 and the corresponding values of n=10,20,40 to cover the interval [0,5]. We label the resulting list of values by A1,A2,A3.

> A1:=Eulermethod(f,0,3,0.5,10);

> A2:=Eulermethod(f,0,3,0.25,20);

> A3:=Eulermethod(f,0,3,0.125,40);

>

Remember, the three lists of pairs correspond to approximate solutions y(x) on [0,5] provided by the Euler method with smaller and smaller step size. Recall that we considered the same IVP : , in Sections 1 and 2. The equation was labeled ODE2 and the numeric solution for the IVP provided by Maple was labeled " soln ". We obtain the plot of a numeric solution with " odeplot " command.

> plo1:=odeplot(soln,[x,y(x)],0..5,color=magenta,scaling=constrained):

> plo2:=pointplot(A1,color=blue,scaling=constrained):

> plo3:=pointplot(A2,color=green,scaling=constrained):

> plo4:=pointplot(A3,color=black,scaling=constrained):

> display([plo1,plo2,plo3,plo4]);

>

As you see, our Euler solutions are getting closer and closer to Maple's solution.

>

Homework Problems

Do your homework in the same Maple session during which you reexecuted the commands in the above worksheet. Do not put the "restart" command in your homework worksheet. If you do, you will have to reload " with(plots) " and " with(DEtools) " as well as copy, paste and reexecute the definition of Eulermethod procedure.

>

Problem 1. Consider the following ODE

(a) Find the general solution to the equation.

(b) Plot the slope field corresponding to the equation for x and y between -2 and 2.

, ,

, .

(d) Plot the two solutions and the slope field in one coordinate system for x between -2 and 2 using the DEplot command.

>

Problem 2. The growth of a certain animal population is governed by the equation

where denotes the number of animals at time t, t is measured in months. There are 200 animals initially.

(a) Find .

(b) Graph in a large enough interval to see the longterm behavior of the population.

Hint : Remember to define the solution to your ODE as an expression or a function, as in Example 3, before attempting to process it.

>

Problem 3. Consider the IVP

, .

(a) Solve the IVP numerically using Maple's numerical solver.

(b) Plot the solution in the interval [0,4] using the "odeplot" command.

(d) Display your Euler approximations and the Maple's solution in one coordinate system.

>

Last modified August 1999.

>