INTRODUCTION TO MAPLE
II. ELEMENTARY CALCULUS AND
DIFFERENTIAL EQUATIONS

Richard M. Felder
Department of Chemical Engineering
North Carolina State University

This tutorial outlines how to use MAPLE to differentiate functions, integrate functions analytically and numerically, and solve systems of first-order ordinary differential equations. To work through it, you should have a good working knowledge of MAPLE command syntax, arithmetic operations, defining expressions and functions (and what the difference is between them), and plotting. If you don't, you are advised to work through the first tutorial in this series and then come back to this one.

TUTORIAL

Open MAPLE and create a worksheet file

If you're running on EOS, bring up the application menu (point to an open space on the desktop and click the middle mouse button) and choose Math and Statistics - Maple(4). (Alternatively, go into the directory where you wish to work and at the command prompt in the Xterm window, type add maple and then xmaple &)

A MAPLE window should come up, with an untitled worksheet superimposed over it.

• Choose Save as under the File menu, and name the file maple2.mws

• Type the indicated commands (shown in boldface type next to the MAPLE prompt [>) followed by the "Return" or "Enter" key.

Tip:

Once you've done something in a program step (like redefine a variable), it stays done unless you change it again. Even if you click on an earlier step and re-execute it, the effect of the change remains.

If a program does not seem to be working, before trying anything else re-execute all steps from the beginning by clicking on the Edit menu and selecting Execute - Worksheet. (You can also select a set of commands and re-execute them by selecting Execute - Selection.)

When completing this tutorial, type the boldface commands preceded by the MAPLE prompt [> (don't retype the prompt), and hit the "Return" or "Enter" key after typing the semicolons or colons that conclude each command. Following each command that ends with a semicolon, we will give the MAPLE output in italics. (If a command ends with a colon the output is suppressed.)

DIFFERENTIATION

Differentiating an expression (as opposed to a function).

Once an expression for f(x) has been defined, use the   diff(f,x)   command to obtain an expression for the first derivative (df/dx), the   diff(f,x,x)   command to obtain the second derivative (d2f/dx2), etc.

In the set of commands that follows, we define an expression for a variable, f(x); review how to evaluate it for x=1, determine expressions for its first and second derivatives, and evaluate those expressions with four significant figures for the same value of x.

[> restart;

[> f   :=   x^3 + 5*exp(2*x);

f   :=   x3 + 5e2x

[> subs(x=1,f);

1 + 5e(2)

[> evalf(subs(x=1,f), 4);

1.677

[> df   :=   diff(f,x);

df   :=   3x2  10e(2x)

[> df1   :=   evalf(subs(x=1,df),4);

df1   :=   1.647

[> d2f   :=   diff(f,x,x);

d2f   :=   6x + 20e(2x)

[> evalf(subs(x=1,d2f),4);

8.706

• Save

Differentiating a function.

Even though an expression for f(x) and a function f(x) may look the same to you, they are treated in totally different ways by MAPLE. In this subsection, we define a function g(x), determine its first and second derivatives using the D(g) and D(D(g)) commands, and evaluate the derivatives at x=1. Remember that the diff command is used to differentiate expressions and the D command is used to differentiate functions.

[> restart;

[> g   :=   x -> x^3 + 5*exp(2*x);

g   :=   x -> x3 + 5e(2x)

[> dg   :=   D(g);

dg   :=   x -> 3x2  10e(2x)

[> evalf(dg(1), 4);

1.647

[> d2g   :=   D(D(g));

d2g   :=   x -> 6x + 20e(2x)

[> evalf(d2g(1),4);

8.706

• Save

INTEGRATION

Indefinite integration of an expression

The command   int(f,x)   is used for both indefinite and definite integration of an expression (not a function) for f(x).

The next commands define an expression for f(x), integrate it (the expression for the integral does not include a constant of integration), evaluate the integral [F(x)] at a specified value of x, and check the integration by differentiating F to obtain f again. Note that MAPLE is case-specific, so that it considers f and F to be two different variables.

[> restart;

[> f   :=   3*x^2  10*exp(2*x);

f   :=   3x2  10e(2x)

[> F   :=   int(f,x);

F   :=   x3 + 5e(2x)

[> evalf(subs(x=1,F), 4);

1.677

[> diff(F, x);

3x2  10e(2x)

We next convert the expression for F(x) into a function using the unapply command, evaluate the function at x=1, and differentiate F(x) with the D(f) command to obtain the original expression for f(x) as a function.

[> F   :=   unapply(F,x);

F   :=   x > x3 + 5e(2x)

[> evalf(F(1),4);

1.677

[> D(F);

x > 3x2  10e(2x)

• Save

Definite integration

The int command may also be used to calculate a definite integral. Instead of int(f,x), the command becomes int(f,x=a..b), where a and b are the lower and upper limits of the integral.

[> restart;

[>f   :=   3*x^2  10*exp(2*x);

f   :=   3x2  10e(2x)

[> A   :=   int(f, x=0..1);

A   :=   4 + 5e(2)

[> evalf(A,4);

 3.324

• Save

Numerical integration

Some functions cannot be integrated analytically, but the definite integrals of such functions still have meaning (the area under a plot of f vs. x between the limits of integration) and MAPLE can determine them. For example, the function   f(x) = exp(  x3)   has no analytical integral (convince yourself), but a plot of f vs. x> is a well-behaved curve that looks like this:

(The apparent gap near the f axis is not real--the curve goes smoothly to f=1 at x=0.)

Let us integrate this function from x=0 to x=2.

[> restart;

[>f   :=   exp(  x^3);

f   :=   e(  x3)

[> F   :=   int(f, x=0..2);

MAPLE returned only the definite integral, indicating that it could not find an analytical expression. However, it has the desired value.

[> evalf(F,4);

.8394

• Save

SOLVING ORDINARY DIFFERENTIAL EQUATIONS

The equations to be solved have the following general form:

where

t is the independent variable

y1(t), y2(t),..., yn(t) are the dependent variables (or state variables)

y1i is the initial condition for the first state variable (the value when t=to)

f1, f2,...,fn are known functions of t and the y's (state functions).

MAPLE can either solve the differential equations analytically (if analytical solutions exist) or numerically. The solutions are expressions for y1(t), y2(t),...,yn(t) (if the equations are solved analytically) or plots or tables of y1, y2,..., yn at values of t from t= to (which is usually 0) to some final value.

Example

A pair of chemical reactions, A->B and B->C, take place in a batch reactor, starting with pure A at a concentration CA0 = 1.00 mol/liter. The following equations describe how the concentrations CA(t), ,CB(t), and CC(t) vary with time(sec). The notation CA' will be used to represent the derivative dCA/dt.

CA' = -0.1CA ,   CA(0) = 1

CB' = 0.1CA - 0.2CB ,   CB(0) = 0

CC' = 0.2CB ,   CC(0) = 0

These equations have the form of the set that began this document. (The independent variable, t, does not appear explicitly in the state functions on the right-hand sides, but nothing says it has to.)

Initially the reactor contains pure A at a concentration of 1.00 mol/L. As time increases, A is steadily consumed in the first reaction . B forms in the first reaction and its concentration initially increases, but it is also steadily consumed in the second reaction and so its concentration eventually starts to decrease and approaches zero at large times. C forms in the second reaction, slowly at first since it cannot appear until some B has been formed, and eventually it is all that is left in the reactor. You would expect plots of the three concentrations vs. time to appear as follows:

The curves of Ca vs. t and Cb vs. t are not identical at large times but they are too close together to be distinguished on the scale of this plot.

The goal is to write a MAPLE program to solve the equations and to calculate and plot Ca, Cb, and Cc at times from t=0 to t=40 s. First solve the equations analytically, then numerically.

[> with(plots): (return key)

This command makes a number of plotting routines available to you. (You could see what they are by using a semicolon instead of a colon at the end of the command.) We will use the plot, textplot, and display routines in the package at different places in the tutorial.

Enter the differential equations to be solved

When entering the next command, don't worry about the error messages that appear after you hit the return key without typing a semicolon or colon at the end of the first three lines; just keep typing. MAPLE will treat the second, third, and fourth lines as continuations of the statement.

[> deqs   :=   (return key)

> diff(Ca(t),t) = -(1/10)*Ca(t) , (return key)

> diff(Cb(t),t) = (1/10)*Ca(t) - (2/10)*Cb(t), (return key)

> diff(Cc(t),t) = (2/10)*Cb(t) ; (return key)

You could have entered the MAPLE text on one line, but it is easier to read and modify the statement if you enter each equation on a separate line.

• Save.

Create the set of initial conditions and a list of dependent variable names

[> inits   :=   Ca(0)=1 , Cb(0) = 0 , Cc(0) = 0;

inits   :=   Ca(0)=1 , Cb(0) = 0 , Cc(0) = 0

[> fcns   :=   [Ca(t), Cb(t), Cc(t)];

fcns   :=   [ Ca(t), Cb(t), Cc(t) ]

• Save.

If you wished to solve the equations numerically, you would skip the next few sections. You may do so now or continue with the analytical solutions.

Solve the equations analytically (symbolically)

[> dsolve({deqs, inits}, fcns);

The command dsolve solves the set of differential equations defined as deqs (you can choose any label) with initial conditions defined as inits to determine expressions for the dependent variables defined as fcns. It has the capability of solving equations of higher order as well as first-order equations. Also, by adding optional arguments in the calling statement, you can get solutions using Laplace transforms (method=laplace), series solutions of any desired order (type=series), or numerical solutions using a variety of algorithms (type=numeric).

• If an analytical solution cannot be found for any reason, MAPLE does not display anything but simply goes to the next prompt. You can anticipate this happening if any of the equations are nonlinear. If it happens, first make sure all of your commands are correct. If they are, you will have to solve the equations numerically. Scroll down to find the procedure for doing so.

• Unless you are obtaining a numerical solution, you must use exact arithmetic in writing the equations that you will use dsolve to solve. The equations may include integers [5], functions of integers [exp(5)], and rational fractions (1/10), but not decimals. If you had put 0.1 in the above commands instead of (1/10), dsolve would not have worked.

• For more details about the command, you can type ?dsolve . (Don't do it now.)

Create functions for Ca(t), Cb(t), and Cc(t)

Even though dsolve seems to have solved the equations and printed the solutions, MAPLE does not "know" what the expressions for Ca(t), Cb(t), and Cc(t) are. To do anything with the solutions, like plotting them, you need to first define them as functions, which is what we do next.

For each of the following commands, you can copy and paste from the dsolve output instead of retyping the expressions. (Use the "Copy" command in the Edit menu, not "Copy as Maple text.")

[> Ca   :=   t -> exp( -1/10*t) ;

[> Cb   :=   t -> exp( -1/10*t) - exp( -1/5*t) ;

[> Cc   :=   t -> 1 - 2*exp( -1/10*t) + exp( -1/5*t);

Generate a plot of the concentrations vs. time from t=0 to t=40 s.

The command that follows is a single command, which we break into different lines for ease of viewing. The parenthetical explanations not in boldface should not be typed.

[> P   :=   plot( [Ca(t), Cb(t), Cc(t)], t=0..40,

> title = `CONCENTRATION (MOL/L) VS. TIME(SEC)`, (use forward quotes)

> titlefont = [HELVETICA, BOLD, 14], (font, font style, font size in pts)

> labels = [`t`, `C`], (axis labels - 6 characters max.)

> style=[point, point, point], (point as opposed to line)

> symbol=[circle, box, diamond],

> color=[red,green,blue] ):

[> display ([P]);

A plot like the one shown at the beginning of this tutorial should appear.

• Save

Label the curves on the plot

First add a paragraph before the "display" command.
• Click on the long vertical square bracket ( [ ) to the left of the "display" command and the chart that follows it. The bracket should be highlighted.
• Click on the Insert menu and select Paragraph Before
• Click in the space that opens, then click on the Insert menu and select Maple Input. A > prompt should appear.

Next prepare labels for the three curves on the chart. With the cursor next to the prompt (>) on the line just created, type on a single line

labela   :=   textplot( [.8, .5, `Ca`] ):
labelb   :=   textplot( [4, .26, `Cb`]):
labelc   :=   textplot( [20, .8, `Cc`] ):

followed by the "Return" key.

The variable "labela" consists of a label (Ca) and the coordinates (abscissa = 0.8, ordinate = 0.5) of the desired location of the label on the plot to be generated by the subsequent display command.   labelb and labelc have similar meanings. It generally takes a few tries to get the text exactly where you want it.

Finally, modify the "display" command to insert the labels. Make the command read

[> display([P, labela, labelb, labelc]);

Then hit the "Return" key to generate the modified chart.

• Save

We could at this point generate a table of values of Ca, Cb, and Cc at specified times, but we will defer that until we have generated numerical solutions.

The linear first-order differential equations in this example have an analytical solution, which MAPLE found. Often analytical solutions to differential equations do not exist and numerical equations must be obtained. In this section we show how to do so.

First, turn Ca, Cb, and Cc from functions back into variable names.

[> Ca := 'Ca':   Cb := 'Cb':   Cc := 'Cc':

If you had not generated analytical solutions first but simply wanted to get numerical solutions, you would have jumped right from the fcns :=  ... statement to the statements in the next section.

Generate numerical solutions of the equations and plot them

To solve differential equations numerically, you must first add MAPLE's differential equations package.

[> with(DEtools):

If you typed a semicolon instead of a colon, you would get a list of the different routines in this package. The only one we will use is DEplot.

Next, create and store separate plots of the dependent variables vs. the independent variable (but don't display them yet).

[> Caplot   :=   DEplot ( [deqs], fcns, t=0..40, [[inits]], scene=[t,Ca], linecolor=red, stepsize=0.2, arrows=none):

[> Cbplot   :=   DEplot ( [deqs], fcns, t=0..40, [[inits]], scene=[t,Cb], linecolor=green, stepsize=0.2, arrows=none):

[> Ccplot   :=   DEplot ( [deqs], fcns, t=0..40, [[inits]], scene=[t,Cc], linecolor=blue, stepsize=0.2, arrows=none):

If you get an error message after any of these commands, go back and re-execute the deqs :=  , inits :=  , and fcns :=   commands and try again.

• Save

The first of the DEplot commands leads to the numerical solution of the equations in deqs with initial conditions in inits. The solution will be values of the variables in fcns (Ca, Cb, and Cc) at intervals of 0.2 time units (stepsize) from t=0 to t=40. After the command is executed, Caplot will contain the coordinates of a plot of Ca vs. t (the variables specified in the scene command) that will appear as a solid red curve when the display command is used to show the graph. The "arrows=none" specification suppresses the generation of a direction field. (If you don't know what that means, don't worry about it.)

The DEplot command uses a fourth-order Runge-Kutta algorithm (other options exist if you specify them with additional arguments) to generate plots of solutions of differential equations. If you ended these three statements with semicolons, MAPLE would have displayed the coordinates of the plots at 0.2 second intervals along with other information about the plot, but not the plot itself. If you just entered the DEplot commands (rather than Caplot   :=   DEplot...) followed by semicolons, the three plots would have been displayed on separate graphs.

By adding options to the DEplot command you can do a number of additional things, like specifying display ranges for the dependent variables, stopping the integration if the dependent variables go outside the specified ranges, changing solution algorithms, generating direction (slope) fields, and integrating higher-order differential equations. To see the options, click on "DEplot" in any of the commands where it appears and then select the HELP menu, or type [> ?DEplot . (Don't do it now.)

Next, if you have not already done so in this tutorial, prepare labels for the three curves on the chart. Click below if you previously created labela, labelb, and labelc.

Label the curves on the plot

With the cursor next to the prompt (>) on the line just created, type on a single line

labela   :=   textplot( [.8, .5, `Ca`] ):
labelb   :=   textplot( [4, .26, `Cb`]):
labelc   :=   textplot( [20, .8, `Cc`] ):

followed by the "Return" key.

The variable "labela" consists of a label (Ca) and the coordinates (abscissa = 0.8, ordinate = 0.5) of the desired location of the label on the plot to be generated by the subsequent display command.   labelb and labelc have similar meanings. It generally takes a few tries to get the text exactly where you want it. We are finally ready to display the three generated plots on a single chart.

Display the graph

[> display ( [Caplot, Cbplot, Ccplot, labela, labelb, labelc] );

The desired chart should appear. If the plots appeared as a set of segments with discontinuous slopes instead of smooth curves, we would re-execute the DEplot commands with a smaller step size.

• Save

Generate and store numerical solutions of the three equations

[> soln   :=   dsolve({deqs, inits}, fcns, numeric);

soln   :=   proc(rkf45_x) ... end

At this point the three equations have been solved using the 4th-order Runge-Kutta algorithm to generate a numerical function soln(t) containing all three concentrations as functions of time.

Display the concentrations at a specified time (t=5), each with four significant figures.

[> evalf(soln(5), 4);

[ t = 5.,   Ca(t) = .6065, Cb(t) = .2387, Cc(t) = .1548) ]

Display a table of concentrations at times 0, 2, 4,..., 20 with four significant figures. (The next command illustrates MAPLE's powerful programming capability, which we have not discussed in these tutorials.)

[> for k from 0 to 10 do evalf(soln(2*k), 4); od;

[ t = 0,   Ca(t) = 1, Cb(t) = 0, Cc(t) = 0]

[ t = 2,   Ca(t) =.8187, Cb(t) =0.1484, Cc(t) = .03286 ]

[ t = 4,   Ca(t) =.6703, Cb(t) =0.2210, Cc(t) = .1087]

[ t = 20,   Ca(t) =.1353, Cb(t) =0.1170, Cc(t) = .7476]

This repetition statement (or for loop) instructs MAPLE to execute the evalf command repeatedly, substituting k = 0, 1, 2,..., 10 wherever k appears. When k=0, 2*k=0 and MAPLE executes the command evalf(soln(0),4) to evaluate the three concentrations at t=0; when k=1, 2*k=2 and MAPLE executes evalf(soln(2),4), etc. If several commands are placed between "do" and "od;", MAPLE will execute them in sequence each time it goes through the loop.

The tabulated results indicate that Cb goes through a maximum somewhere between 6 and 8 seconds (verify this statement). To get a better estimate of the maximum you could evaluate soln at closely spaced times between t=6 and t=8. (Don't do it now.)

Finally, clean up the worksheet by deleting the next-to-last paragraph. Click on the bracket ( [ ) to the left of the evalf(soln(5),4); command. Choose Delete Paragraph from the Edit menu. The command and the MAPLE response should be gone.

• Save

Insert text in the worksheet to prepare a report (if one is required as an assignment)

• Insert a new paragraph and title at the top of the worksheet

• Click on the bracket ( [ ) to the left of the first command line ("restart;"). The bracket should be highlighted. Click on the Insert menu and select Paragraph Before. A line should open up before the command line.

• On the new line, type

MAPLE TUTORIAL II. CALCULUS AND DIFFERENTIAL EQUATIONS

Don't hit "Return" yet.

• Click on the triangle to the right of the style droplist (the horizontal window with the word "Normal" in it at the left of the MAPLE context bar) to bring up a list of predefined styles. Then click on Title. Your title should now be centered, underlined, and in large type.

• Put your name (or names and group number) as the author(s) of the report

• With the cursor at the end of the title, hit the "Return" key. Type your name(s) or, if applicable, group number, colon, and last names separated by commas. Don't hit "Return" yet.

• Change the style of this line from title to Heading 1 (click on the style droplist button as before). Then use the font size droplist to change the size from 18 point to 14 point. Then click on the "Centered Text" button (second from the right on the context bar). What you typed should now be boldface and centered, in larger type than normal size but smaller than in the title style.

• Hit the "Return" key. You should now be back in Normal style.

• By clicking and dragging the mouse, highlight the input and output regions from the first restart; command through the evalf(d2g(1), 4); command and its output region (8.706) and the Save command immediately following. Click on the Format menu and select Indent. Click to the right of the little square that appears and type DIFFERENTIATION.

• Highlight the region from the first restart; command through the evalf(subs(x=1,d2f), 4); command line and its output region (8.706). Click on the Format menu and select Indent, click to the right of the little square and type Differentiating expressions.

• Highlight the region from the second restart command through the evalf(d2g(1), 4); command and its output region (8.706). Click on the Format menu and select Indent, click to the right of the little square and type Differentiating functions. Save.

• Define sections for the remainder of the worksheet

• Using the same sequence of steps (highlight, indent, type heading), define a section containing the commands from the third restart; command to the evalf(F,4); command and its output region (.8930). Give this section the heading Integration.

• Define a subsection from the third restart; to the D(F); command and its output region (x --> 3x2...). Call it Indefinite integration .

• Define a subsection from the next restart; to the evalf(A,4) command and its output region (-3.324). Call it Definite integration .

• Define a subsection from the next restart; to the evalf(F,4); and its output region (.8930). Call it Numerical integration.

• Save.

• Define a main section from the final restart; to the end of the program, and give it the heading ORDINARY DIFFERENTIAL EQUATIONS. Save.

• Define a subsection from the final restart; to the fcns   :=   [Ca(t), Cb(t), Cc(t)]; command and its output region. Call it Define equations, initial conditions, and list of dependent variable names. Save.

• Define a subsection from the dsolve([deqs,inits],fcns); command to the Cc   :=   t > 1-2*exp(...); command and its output region. Call it Solve equations analytically and create functions for the individual solutions. Save.

• Define a subsection containing the commands from with(plots): command line to the plot that follows the display[P,...]; command, giving it the heading Plot the solutions and label the curves. (You may have to play around with the mouse to highlight the chart but not the command that follows it.) Save.

• Define a subsection containing the commands from with(DEtools): through the chart that follows display...; Give it the heading Generate and plot numerical solutions. Save.

• Define a subsection containing the commands and output regions from soln   :=  ...; to Cc(t) = .7476. Give it the heading Solve equations numerically and generate a table of the solutions. Save.