prepared by Mark E. Johnson, maejohns@math.princeton.edu

MAE 305

Mathematics in Engineering I


Mathematica Tutorial

An outline of this document follows:
Introduction
Basics
Getting Started
Easy Examples
Symbolic Mathematics
Solving differential equations
More to Come
Appendix : Newton-Raphson Code
Power Series Solutions to Differential Equations

Introduction

One of the primary objectives of the MAE305 course is to provide you with exposure to a bag of tricks used in solving differential equations (DEs). Included in this bag of tricks are a number of numerical and symbolic routines which can help you solve DEs on a computer.

The purpose of this Mathematica tutorial and the similar MATLAB tutorial is to provide you with a starting point for the use of these powerful and well-developed programs -- again we hope to provide you with some exposure to alternative methods of solving DEs. In no way are these tutorials intended to be a complete reference to the powers of this software, nor are these two pieces of software the only such software available out there. One that comes to mind is Maple -- equally deserving of a course introduction, but we thought we'd stop after two tutorials. If you want to learn more about Maple, feel free to ask and we can point you in the right direction -- it is a really nice piece of software, too!

To get things started, we are going to jump ahead to the MATLAB tutorial first, then return here to review Mathematica. So if you want, please follow this link to the MATLAB page.

If you have just returned from the Matlab page, then you will notice some redundant information below concerning setting your DISPLAY and all that jazz. You should be able to read through it and skip over the things which you have already done. I have tried to make each of the tutorials self-contained so if you ever want to come back to play with just one or the other, you will be able to just look on the same page to remind yourself how to set everything up...

Now, what are we going to do with Mathematica? Specifically, in this tutorial you will learn to do the following tasks in Mathematica:

This tutorial is very far from illustrating the full capabilities of the Mathematica software -- the intention is to be able to provide you with something to get things started. Once you understand some of the basics, you can come back to this page and look at a more complete tutorial available from Los Alamos National Laboratory ; you also might want to check out the home web page of Wolfram Research, the creator of Mathematica.

For information on using Mathematica at Princeton, follow this link.

Basics

There are a few useful things that you should remember about Mathematica in general. The first one is that all Mathematica commands start with capital letters (this will become obvious in the examples). Next, the way to get on-line help from Mathematica is by using ?CommandName at the prompt. For example to find information on integrating functions we use,
In[1]:= ?Integrate
Integrate[f,x] gives the indefinite integral of f with respect to x.
   Integrate[f,{x,xmin,xmax}] gives the definite integral.
   Integrate[f,{x,xmin,xmax},{y,ymin,ymax}] gives a multiple integral.
You should also remember that the percent symbol (%) used in command arguments will be referring to the expression in the last command line or, if used with a number (e.g. %18), will refer to the expression of that input line number (e.g. Integrate[%18]). Lastly, arguments to the Mathematica command lines are put in square brackets and lists or arrays are put in curly brackets. Oh, and to quit Mathematica use ``Quit''.

Getting Started

Before continuing with the rest of this tutorial, make sure that you have (a) successfully logged in (in this example, replace ``kerouac'', which is my machine's name, with the name of the machine you are actually working on) and (b) created a window (besides the console window) on your screen. To create a new window, go to the ``Tools'' bar in the toolbox and find the ``Shell'' option, OR hold down the right mouse button while the pointer is inside your console widow, drag the mouse to the ``clone'' menu and choose a size for your new shell ; release the mouse button. Now, go to the console window and type,
maejohns@kerouac:~ : xhost + phoenix
phoenix being added to access control list
You should have typed ``xhost + phoenix'' (everything else in the first line
maejohns@kerouac:~ : 
is my prompt). The xhost command gives permission to phoenix to open windows on your machine. Next, go to your other window and type
maejohns@kerouac:~ : rlogin phoenix
Once you are logged in, make a directory called ``MAE305'' in which you will store class information. You can do this by typing the command
mej@phoenix.Princeton.EDU:~ : mkdir MAE305
NOTE: You need only do this once. If you have already created this directory, you will get a ``mkdir: MAE305: File exists'' error message. Now, chang into the directory by issuing the cd command :
mej@phoenix.Princeton.EDU:~ : cd MAE305
Finally, you need to set the DISPLAY environment variable to the name of the machine at which you are working --- this tells phoenix where to open new windows. Don't forget the :0. Let's also set the PRINTER variable.
mej@phoenix.Princeton.EDU:~/MAE305 : setenv DISPLAY YourMachineName:0
mej@phoenix.Princeton.EDU:~/MAE305 : setenv PRINTER PrinterName
Now, some of you (particularly those who have taken CivE 201 in the past) might have trouble with the previous commands. If you have problems, then try first issuing the command ``csh'', then trying the setenv commands. If this works, then you will have to type this each time you plan on setting the DISPLAY variable. If you are still having problems, then please ask me for further help!

The Indigo2's in EQuad Room E423 have YourMachineName listed on a piece of paper on the monitor. You can also use the UNIX command hostname to find out the name of your machine. PrinterName should be replaced with the name of the printer in the computer room in which you are working (or where you want hard-copies to be printed). The printer in EQuad Room E423 is currently named gutenberg. Now, before we proceed with our Mathematica demonstration, copy a file from my home directory to put into your class directory.

mej@phoenix.Princeton.EDU:~/MAE305 : cp ~mej/MAE305/data .
Don't forget to include the period at the end of this command! What this command does is copies the file data from the directoy called ~mej/MAE305 and places it in the current directory (which should be ~yourusername/MAE305). The files name stays the same, namely data. If you wanted to call it something else, then replace the period with a new filename.

We will get back to where this file came from in a little while. The UNIX command ls can be used to list the files in the current directory ; by issuing this command, you should see data listed on the screen.

Now we are ready to start Mathematica by typing the command math.

mej@phoenix.Princeton.EDU:~/MAE305 : math
If you happen to get the following message
mej@phoenix.Princeton.EDU:~/MAE305 : math
math: Command not found.
then let me know so we can make some repairs to some of your initialization files (in particular, to your PATH, which is probably set in your .cshrc or your .login files). For the time being, you can still proceed by typing
mej@phoenix.Princeton.EDU:~/MAE305 : /usr/princeton/bin/math
I would like to point out that there is also another way to start Mathematica from these computers, that is, by typing mathematica instead of math at the prompt. This command will open a window which provides a nice interface for inputting and retrieving data. We, however, will not be using this interface for this session ; if you are interested in learning how to use it on your own, then I highly suggest to do so! Feel free to come by my office during office hours and we can experiment with it together.

Easy Examples

Now that you are in Mathematica, I would like to jump immediately to an example of using it to solve a differential equation

y' - (4x)/(x*x+1) * y = Exp(x) * (1 + x*x) ^ 3

mej@phoenix.Princeton.EDU:~/MAE305 : math
Mathematica 2.2 for SPARC
Copyright 1988-94 Wolfram Research, Inc.
 -- Motif graphics initialized -- 

In[1]:= sol=DSolve[y'[x]-4x/(x^2+1) y[x]==(1+x^2)^3 Exp[x],y[x],x]

                   x       2 2             2          2 2
Out[1]= {{y[x] -> E  (1 + x )  (3 - 2 x + x ) + (1 + x )  C[1]}}
Nice! Now let's slow down and ease into Mathematica with some more elementary demonstrations. First, let us create a vector of numbers from 0 to 2 with a step of 0.1.
In[1]:= x = Range[0.0,2.0,0.1]

Out[1]= {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 
 
>    1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.}
Notice how the result was displayed (echoed) by Mathematica. If you don't want this to happen, put a semicolon ``;'' at the end of the command line,
In[1]:= x = Range[0.0,2.0,0.1];
Now, let's say that we want to look at the square root of x in this range,
In[2]:= f = Sqrt[x];
Now, in order to plot this (f versus x) we have to create a list that Mathematica will understand (i.e. pairs of numbers). First, do the following,
In[3]:= p = {x,f}

Out[3]= {{0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 
 
>     1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.}, 
 
>    {0., 0.316228, 0.447214, 0.547723, 0.632456, 0.707107, 0.774597, 
 
>     0.83666, 0.894427, 0.948683, 1., 1.04881, 1.09545, 1.14018, 1.18322, 
 
>     1.22474, 1.26491, 1.30384, 1.34164, 1.3784, 1.41421}}

You should notice that now we just have two long lists of numbers, not pairs of numbers specifying a point. We have to take the Transpose of this list,
In[4]:= p = Transpose[p]

Out[4]= {{0., 0.}, {0.1, 0.316228}, {0.2, 0.447214}, {0.3, 0.547723}, 
 
>    {0.4, 0.632456}, {0.5, 0.707107}, {0.6, 0.774597}, {0.7, 0.83666}, 
 
>    {0.8, 0.894427}, {0.9, 0.948683}, {1., 1.}, {1.1, 1.04881}, 
 
>    {1.2, 1.09545}, {1.3, 1.14018}, {1.4, 1.18322}, {1.5, 1.22474}, 
 
>    {1.6, 1.26491}, {1.7, 1.30384}, {1.8, 1.34164}, {1.9, 1.3784}, 
 
>    {2., 1.41421}}

Now we can plot this.
In[5]:= ListPlot[p,PlotJoined->True]

Out[5]= -Graphics-

``PlotJoined'' is one of the Mathematica plotting options (or attributes) that we change to make our plots look nice. A list of other options can be found in Appendix B of this handout, and a couple more are demonstrated in this tutorial. Notice that we could have accomplished the same thing by just plotting the function in that range,
In[6]:= Plot[Sqrt[x],{x,0,2},AxesLabel->{"x","Sqrt(x)"}]

Out[6]= -Graphics-
If we want to make a hardcopy of this plot in the file ``outfile'', do the following:
In[7]:= Display["outfile",%6]

Out[7]= -Graphics-

The Display command takes the graphics that was generated in input line number 6, from above, and write it to a file called ``outfile''. However, this file is of no use to us yet, as it is not in PostScript format, and thus would not be understood by the printer. In order to convert this file to a PostScript file, type the following:

In[8]:= !psfix outfile > outfile.ps
Now we have a PostScript file which we can print out by using the ``lpr'' command (``lpr outfile.ps''). Note : the command psfix is actually a UNIX command and can be executed from your main shell window (when you are not running Mathematica), but we see that we can conveniently execute UNIX commands from within Mathematica by placing an exclamation point (!) in front of it. If you were to leave the exclamation point off, then you should get an error. If I am confusing you, then either ask someone (me?) to explain it to you or just ignore what I am saying altogether and remember to use ``!'' while in Mathematica and to leave it off while outside of Mathematica.

Plotting input files

Let us now practice reading in and plotting data from a file. The data which you should have copied from my directory contains the first few iterates from Newton's method to find the square root of 2. A little C program was written to perform a Newton-Raphson method to find a solution to the equation x*x=2 which should be the square root of 2. The Newton-Raphson method can be written as an iterative scheme :

x(n+1)=x(n)-f(x(n))/f'(x(n)),

where f(x)=x*x-2. The file that you copied from my directory at the beginning of the session, ``data'', contains the results of this iteration from an initial guess of x(0)=10. We want to read this data into Mathematica and plot it. The data is in two-column format: n x(n).

In[1]:= data=ReadList["data",{Number,Number}];
In[2]:= ListPlot[data,PlotJoined->True,PlotRange->{{0,10},{0,10}}]
Now edit a file (in another window with emacs, vi, jot, or whichever editor you are most familiar) called mymath.m with the commands which follow. If you don't know how to make such a file, just copy the one I made by typing
cp ~mej/MAE305/mymath.m ~/MAE305/mymath.m
You can view the contents of a file with the cat command: try it by typing cat mymath.m at your UNIX prompt. The contents should read:
data = ReadList["data",{Number,Number}];
ListPlot[data,PlotJoined->True,PlotRange->{{0,10},{0,10}}]
This is our first Mathematica macro. To load in into your session, type
In[1]:= << mymath.m

Symbolic Mathematics

In this section, I will quickly run through some examples displaying Mathematica's symbolic computation power. We can simplify, expand, or factor algebraic expressions
In[1]:= (x + 2y + 1)(x-2)^2



Out[1]= {4. (1. + 2 y), 3.61 (1.1 + 2 y), 3.24 (1.2 + 2 y), 2.89 (1.3 + 2 y), 
 
>    2.56 (1.4 + 2 y), 2.25 (1.5 + 2 y), 1.96 (1.6 + 2 y), 1.69 (1.7 + 2 y), 
 
>    1.44 (1.8 + 2 y), 1.21 (1.9 + 2 y), 1. (2. + 2 y), 0.81 (2.1 + 2 y), 
 
>    0.64 (2.2 + 2 y), 0.49 (2.3 + 2 y), 0.36 (2.4 + 2 y), 0.25 (2.5 + 2 y), 
 
>    0.16 (2.6 + 2 y), 0.09 (2.7 + 2 y), 0.04 (2.8 + 2 y), 0.01 (2.9 + 2 y), 
 
               -31
>    1.97215 10    (3. + 2 y)}

OH NO!! This is not right! Do you remember above where we used the variable x to represent a list (with the Range command)? Well, Mathematica still thinks that x is a list, when we actually want it to be a variable. How do we fix this? Well, there are two ways to get around this. The first, most obvious one, is to not use the same variable name, but use something like x1 instead.


In[2]:= (x1 + 2y + 1)(x1-2)^2

                 2
Out[2]= (-2 + x1)  (1 + x1 + 2 y)


Or we can reset the value of x by using the Clear command as follows:

Clear[x]

So let's now start over


In[1]= (x + 2y + 1)(x-2)^2

                2
Out[1]= (-2 + x)  (1 + x + 2 y)

In[2]:= Expand[%]

               2    3                    2
Out[2]= 4 - 3 x  + x  + 8 y - 8 x y + 2 x  y

In[3]:= Simplify[%]

                2
Out[3]= (-2 + x)  (1 + x + 2 y)
In the following line, a %2 is used to denote the 2nd output (Out[2], from above) in this session. This is probably different from the corresponding number in your session!

In[4]:= Factor[%2]

                2
Out[4]= (-2 + x)  (1 + x + 2 y)

or find derivatives of functions,
In[5]:= D[%4,x]

                2
Out[5]= (-2 + x)  + 2 (-2 + x) (1 + x + 2 y)

In[6]:= D[%4,y]

                  2
Out[6]= 2 (-2 + x)
We can also use Mathematica to help us play with the Partial Fractions game, which will be important in solving differential equations,
In[7]:= Apart[(x^2+8)/(x^2-5x+6)]

              17       12
Out[7]= 1 + ------ - ------
            -3 + x   -2 + x

or to solve systems of equations,
In[8]:= Solve[{ax+y==0,2x+(1-a)y==1},{x,y}]

               1   (1 - a) ax
Out[8]= {{x -> - + ----------, y -> -ax}}
               2       2
or to perform symbolic integration,

In[9]:= Integrate[(x^2+8)/(x^2-5x+6),x]

Out[9]= x - 12 Log[2 - x] + 17 Log[3 - x]
Now, before moving on to solving some differential equations, I would go back to the discussion of Newton's method above. We can avoid using a C program altogether by defining a recursion relation within Mathematica:
In[1]:= y[n_]:= y[n] = y[n-1] - (y[n-1] y[n-1] - 2.)/(2 y[n-1])

In[2]:= y[0]=10

Out[2]= 10

In[3]:= ytable = Table[y[i],{i,0,10}]

Out[3]= {10, 5.1, 2.74608, 1.73719, 1.44424, 1.41453, 1.41421, 1.41421, 
 
>    1.41421, 1.41421, 1.41421}
Note the decimal point after the 2 in the recursion relation. This tells Mathematica to do floating point arithmetic rather than symbolic. Compare these numbers with those given in the ``data'' file from my directory.

Obviously, this is a short list of Mathematica examples ; more commands can be discovered by exploring other tutorials which are on the Web or by browsing through one of the many Mathematica books in the library (at least one of which is on reserve).

Solving differential equations

Yes! Mathematica can even be used to solve differential equations. Only in the simplest examples is the computer able to solve the equation without a little help from us. damn. Here is one such example (for which I hope you don't need a computer!).

First we need to clear the values of x and y by typing

Clear[x,y]
Now continue as follows:
In[10]:= equation=y'[x]==y[x]

Out[10]= y'[x] == y[x]

In[11]:= DSolve[equation,y[x],x]

                    x
Out[11]= {{y[x] -> E  C[1]}}
For some more examples, I refer you to the usual places : the Web and particularly the book ``Differential Equations with Mathematica'', by M.L. Abell and J.P. Braselton, Academic Press Professional, 1993.

More to come

This tutorial is an attempt to introduce to you some of the basics of the software package Mathematica. Becoming completely familiar with new software takes time and practice (and a lot of reference to books, friends, and help facilities). Enjoy hacking and don't be afraid to send me email if you have any questions (maejohns@math.princeton.edu).

Appendix : Newton-Raphson Code

#include< stdio.h >
main() {
  double next, last ;
  int i ;
  last = 10.0 ;  /* initial guess */
  for (i=0;i<10;i++) {  /* perform 10 iterates */
    next = last - (last*last-2.0)/(2.0*last) ;
    printf("%d %f\n",i,last);
    last=next;
  }
}

Power Series Solutions to Differential Equations

In this section, we wish to outline how to solve the differential equation
y''+y=0
using power series approximations subject to two different initial conditions.

The first set of initial conditions we wish to use are y(0)=1, y'(0)=0. We will solve this equation approximately on the computer, but since you can solve this problem exactly by hand, we have the ability to compare our approximate solution to the exact one.

The first thing we do in Mathematica is define the left hand side (or lhs) by:

lhs = y''[x] + y[x]
There is no need to define the right hand side since it is zero in this case. Now we want to express this left hand side as a power series in x about the point 0 up to 7th order.
ser = Series[lhs, {x,0,7}]
Now we will take this result and substitute the initial conditions (this is accomplished with the "/." notation).
serone = ser /. {y[0]->1, y'[0]->0}
Now that we have done this, we tell Mathematica to collect similar powers of x after we set the left hand side equal to the right hand side and writing it in a logical expression (no need to know what this means):
equations = LogicalExpand[serone==0]
and now we solve the equations
roots = Solve[equations]
Now that all the coefficients have been solved for, we need to plug them into a series approximation for the solution y[x]. First we define the series
sery = Series[y[x],{x,0,9}]
and now we substitute the initial conditions and the solution "roots" obtained above.
solapprox = Normal[sery] /. {y[0]->1, y'[0]->0} /. roots[[1]]
To plot this solution on the interval from 0 to 3, we type
Plot[solapprox,{x,0,3}]
and to plot both this solution and the ACTUAL solution, we do the following
Plot[{solapprox,Cos[x]},{x,0,3}]
The two solutions seem to match pretty well on this interval. Try plotting them over a largeer interval to see what happens:
Plot[{solapprox,Cos[x]},{x,0,10}]

Different Initial Conditions

Now let's try the same problem with different initial conditions. I will simply include all of the same commands below, but changed to account for the initial conditions y(0)=0, y'(0)=1.
lhs = y''[x] + y[x]
ser = Series[lhs, {x,0,7}]
serone = ser /. {y[0]->0, y'[0]->1}
equations = LogicalExpand[serone==0]
roots = Solve[equations]
sery = Series[y[x],{x,0,9}]
solapprox = Normal[sery] /. {y[0]->0, y'[0]->1} /. roots[[1]]
Plot[solapprox,{x,0,3}]
Plot[{solapprox,Sin[x]},{x,0,3}]
Plot[{solapprox,Sin[x]},{x,0,10}]