Introduction

Basics

Getting Started

Easy Examples

Symbolic Mathematics

Solving differential equations

More to Come

Appendix : Newton-Raphson Code

Power Series Solutions to Differential Equations

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*:

- create data vectors and arrays
- read-in data from files
- plot basic functions (sine, cosine, square-roots, etc)
- plot data
- use the
*Symbolic Manipulation*power of*Mathematica*. and apply it to solving ODEs - make printout of your plots (hard-copies)

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

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

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,

mej@phoenix.Princeton.EDU:~ : cd MAE305Finally, you need to set the

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

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

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

mej@phoenix.Princeton.EDU:~/MAE305 : /usr/princeton/bin/mathI would like to point out that there is also another way to start

**
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

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

In[1]:= x = Range[0.0,2.0,0.1];Now, let's say that we want to look at the square root of

In[2]:= f = Sqrt[x];Now, in order to plot this (

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

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

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

**
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

data = ReadList["data",{Number,Number}]; ListPlot[data,PlotJoined->True,PlotRange->{{0,10},{0,10}}]This is our first

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)}

In[2]:= (x1 + 2y + 1)(x1-2)^2 2 Out[2]= (-2 + x1) (1 + x1 + 2 y)Or we can reset the value of

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[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

In[7]:= Apart[(x^2+8)/(x^2-5x+6)] 17 12 Out[7]= 1 + ------ - ------ -3 + x -2 + xor 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 2or 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

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

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

#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

ser = Series[lhs, {x,0,7}]Now we will take this result and substitute the initial conditions (this is accomplished with the "

serone = ser /. {y[0]->1, y'[0]->0}Now that we have done this, we tell Mathematica to collect similar powers of

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

sery = Series[y[x],{x,0,9}]and now we substitute the initial conditions and the solution "

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}]