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:
For information on using Mathematica at Princeton, follow this link.
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''.
maejohns@kerouac:~ : xhost + phoenix phoenix being added to access control listYou 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 phoenixOnce 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 MAE305NOTE: 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 MAE305Finally, 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 PrinterNameNow, 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 : mathIf 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/mathI 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.
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.psNow 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.
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.mYou 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
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).
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.
#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;
}
}
y''+y=0using 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}]
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}]