Introduction to Mathematica

Richard Gass
Department of Physics,
University of Cincinnati,
Cincinnati, OH 45221-0011
Richard.Gass@UC.EDU

*********************************COPYRIGHT**************************

        Copyright, 1997, 2010, all rights reserved, Richard Gass (Department of 
Physics, The University of Cincinnati, Cincinnati, OH 45221-0011).  This 
material may not be reproduced for profit, modified or published in any 
form (this includes electronic redistribution) without the prior written 
permission of the author listed above. 

**********************************************************************

Getting Started

Mathematica is divided into two pieces: the front end and the kernel. The front end provides the user interface and the kernel does the actual  computations. The front end and the kernel communicate with each other using a program called MathLink. The front end and the kernel need not be on the same machine, you can run the front end on a Macintosh or a PC and run the kernel on (for example) a remote UNIX workstation. The window that this text appears in is called a notebook and is the user interface.

The Menus

If you look at the top of the screen you will see a menu bar. There are  over a hundred menu items (many of which have sub-menus). Listing them all would not be helpful. Instead I have tried to give you an  abbreviated list of the more important commands and to indicate how they are organized. 

The File command gives you a fairly standard list of options: open, close, new, save, save as, print, and print selection plus a few others which we will worry about later. The Edit command gives you the standard cut, paste, and copy options as well as the option to save selection in special format. A few more specialized commands are also available under the cell menu. The Cell menu controls the behavior of cells, while the Format menu controls the appearance of the contents of the cells. The Input menu allows you to make buttons, palettes and hyperlinks. It also allows you to get file paths, graphics coordinates and three dimensional viewpoints. The Kernel menu contains commands for controlling the evaluation of expressions and for selecting the kernel that you want to use. The Find menu has commands for finding and replacing strings in a notebook and for making an electronic index of a notebook. The Window menu controls how the windows for the notebooks that you have open are  displayed. The Help menu takes you to the help browser. The on-line help is first rate.

Running Code

To get started let's type some input. After you have entered the input line, hit either Enter or Shift/ Return, not Return. If your keyboard has two Enter keys and no Return key hit the  Enter key which is on the numerical pad to the right of the keyboard. I will always refer to this key as the Enter key and to the Enter key on the main keyboard as the Return or new line key. The Return  key gets you a new line and nothing else. As you move your mouse through the notebook you will see the mouse pointer change from vertical to horizontal . The pointer is vertical if your mouse is in a cell. Cells are delineated by the blue brackets on the right hand side of the notebook. The pointer is horizontal if your mouse is between two cells (or below the last cell in your notebook). When you mouse pointer is horizontal click once with your mouse. You will get a horizontal line and you can now type in some input. Try adding two numbers. Now, let's try a more complicated piece of input. To run the code below put your cursor anywhere in the code below and hit Enter.

Factor[77/4 - (41 x)/4 + (127 x^2)/6 - (395 x^3)/12 + 
(17 x^4)/4 - (127 x^5)/6 + (41 x^6)/3 + 6 x^7]

Notice that the input is in bold! The output is in plain text just like this sentence. The input line has an input number associated with it (In[1]) and the output line has an output number associated with it (Out[1]).  The command, Factor, does just what it says, it factors polynomials. Notice that the word factor in the previous sentence is blue and underlined. Anything in a notebook that is blue and underlined is a hyperlink. Try clicking one it. Notice that the command Factor starts with a capital letter. All Mathematica commands start with capital letters! The polynomial to be factored is enclosed in square brackets. All Mathematica  functions use square brackets to enclose their argument. If I type (and run):

Factor(x^2 + 2 x + 1)

I get a syntax warning message, a suggested change and the command, Factor, is echoed back. In Mathematica, parentheses are used to change the order of operations and curly braces are used for lists. 
If I try:

factor[x^2 + 2 x + 1]

I get the function echoed.  Mathematica  knows about Factor, but assumes that factor is a user defined function. When you define a function it is a good idea to start its name with a lower case letter to avoid conflict with built in Mathematica functions.

 There are several ways to indicate multiplication in  Mathematica . The first is just like Fortran: use an asterisk.

2*3

Or, you can use a space.

2 3

If you leave out the space then of course you get twenty three. The same thing is true with symbols. Thus, x y and x*y and (x)(y) are all x times y but xy is a new symbol. However, both 2 x and 2x give 2 times x.  This is the only exception to the spacing rule.  Because of this you cannot start a variable name with a number.  Mathematica  assumes that if you have a number in front of a variable, you want them multiplied. Space does not matter for addition or subtraction. That is, a+b and a + b both give the same result. Powers are handled the same way as in Fortran 90: you use ^ for a power. Thus 2^(1/2) is the square root of 2. You can also use Sqrt[2].  However, you can not use ** for  powers in Mathematica. In Mathematica a**b is the non-commutative multiplication not ab. 

Almost all Mathematica  commands are spelled out completely, they are not abbreviated. The exceptions are Abs (absolute value), Re (real part), Im (imaginary part), Sqrt (square root), and D (derivative). Wolfram Research is quite proud of this. Most people dislike this feature at first (it would be quicker for instance to be able to type Int rather than Integrate) but Mathematica has about a 1500 commands and abbreviations would make them very hard to remember. 

You do not have to run Mathematica  commands in the order you originally  typed them in. To run one of your commands, use the mouse to put the blinking cursor anywhere on the desired input command and hit Enter. This feature is useful when you are trying to get code to work or are just playing around, but as we shall see, it can get you into trouble.

Mathematica has all the common mathematical constants built in. Mathematica  uses Pi for π, E for e and I for i. As we will see later you can also use π, e and i.

Getting Help

If you do not remember the form of a command you can ask Mathematica . Suppose you forget the syntax for Plot3D, type  ?Plot3D. Typing ??Plot3D gets you even more information but as you can see, in most cases ? tells you all you want to know. The options information that ?? gets you about Plot3D will be discussed in the chapter on graphics. Clicking on the blue >> symbol will take you to the help pages on the function.

?Plot3D

??Plot3D

If you cannot remember the whole command, type just the part you are sure of. For example, if you type ?Plot*, you will get a list all Mathematica functions that start with Plot. To get the  corresponding syntax type ? followed by the function name as in the example above.

?Plot*

If you type  ?Int* you get:

?Int*

The asterisk is a wild card variable which matches anything.
Thus typing  ?*Pol* gives:

?*Pol*

You can also get help by using the Help browser. To look up the command Integrate go to the Help menu and select Documentation Center. When the Documentation Center comes up, type Integrate in the text block and hit return.

Making an Expression a Number

If you ask Mathematica  for the sine of 1 you get Sin[1] not a number. This is because Mathematica tries to avoid losing accuracy by replacing an exact answer by a decimal approximation. You can get a decimal number by using the function N. If you ask for the sine of 1., you will get a decimal number since 1. is itself a decimal.

Sin[1]

N[Sin[1]]

Sin[1.]

Pi

a=Sin[π/3]

N[Pi]

N[a]

N[π]

N[π,10000]

Adding Text and Changing Styles

Anything that you type in is by default Input and will appear in a bold type face. You can make it text by selecting  the cell bracket (that’s the square bracket on the far right of the page) and then going to the Style menu. Select  Cell Style and then select text from the pull down menu. You can help organize your work by breaking it up into sections and subsections. You make things sections and subsections in the same way. First select the cell, then go to the Style menu and pick Cell Style. Use the pull down menu to change the cell to what ever style you like. You can also use the Style menu to change the face, size or color of your type.

This is text.

Basic Mathematica

Defining Functions

Functions with Immediate Evaluation

Suppose we want to define a simple function: say f(x)=(x^2). The proper way to define such a function in Mathematica is

Clear[f]
f[x_] = x^2

The Clear command removes any previous definition of f(x).  Notice that  I have used x_ this makes x a "wild-card" variable so that anything that is in the slot is matched and squared. Wild-card variables must appear only on the left hand side of a function definition. To see what happens if you fail to use a wild-card variable let's define a second function g(x):

Clear[g]
g[x]=x^2

The function f[x] works correctly

f[z+1]

f[2]

but g[x] does not:

g[x]

g[z+1]

g[2]

Functions with Delayed Evaluation

Sometimes you want to define a function that is  evaluated at the time it is called, rather than when it is defined. To define a function with delayed evaluation  use := rather than =. For example,

Clear[h,a,b,c]
h[x_]:=a x^2+b x+c

Notice that there is no output from a delayed function definition. If I call the function then it evaluates:

h[x]

Plotting

Plotting in Two Dimensions

Plotting is easy in Mathematica. You just use the command Plot. The syntax is Plot[f[x],{x,xmin,xmax}]. For example:

Plot[Sin[x^2] ,{x,0,2Pi}]

If you want to plot more than one function on a single plot, then you give Plot a list (Mathematica uses {} to denote a list).

Plot[{Sin[x^2],Cos[x^2]},{x,0,Pi}]

Of course all plots should be labeled. To label the axes you use the option AxesLabel->{"x axis label","y axis label"}. To label the plot you use the option PlotLabel->"text string". You must enclose the text string in double quotes or you will get the current value of the text rather than the text. All Mathematica options are set using the syntax Option->Value.

Plot[{Sin[x^2],Cos[x^2]},{x,0,Pi},PlotLabel->"sin(x^2) and cos(x^2)",
AxesLabel->{"x","y"}]

There are a number of options for plot. To see what they are type ??Plot

To learn how Plot works try out the function “showPlotPoints”. You do not need to understand how the code works to use this function. The calling syntax is the same as the syntax for plot.

showPlotPoints[f_,{x_,xmin_,xmax_},options___]:=Module[{},{plot,points}=Reap[Plot[f,{x,xmax,xmin},PlotStyle->Black,EvaluationMonitor:>Sow[{x,f}],options]];
Show[plot,ListPlot[points]]]

showPlotPoints[Exp[-(x^2)/3],{x,-3,3}]

The blue points show the points at which the function being plotted was evaluated.

Plotting in Three Dimensions

To plot a function of two variables you use Plot3D. The syntax is Plot3D[f[x,y],{x,xmin,xmax},{y,ymin,ymax}].

Plot3D[Sin[x^2+y^2],{x,-Pi/2,Pi/2},
{y,-Pi/2,Pi/2},PlotPoints->90,Mesh->False]

Like Plot, Plot3D takes many different options. To see the list of options type ??Plot3D.

Additional Plotting Functions

Mathematica  has a number of built in plotting functions. We can ask Mathematica for a list of these functions by typing ?Plot*.

?*Plot*

Not all of these are plotting functions. You can use ?"name" to find out what each of these commands does. There are still other plotting routines that must be loaded as packages. A complete list o packages is given in the reference guide here. To load a package, type Needs["package`"]. To load the package ErrorBarPlots  type:

Needs["ErrorBarPlots`"]

To see what new plotting functions have been loaded type ?*Plot* again. Notice that the new functions are listed separately.

?*Plot*

Exercise Set One

Problem 1.1		Express sin(Pi/7), E^Pi and 7/11 as decimal numbers. Compute the natural log of pi to 1000 places. Compute 400! . 
Note that Mathematica uses Log for the natural log. Log[b,x] gives the log of x to the base b.

Problem 1.2		Define a function to compute the sum of two numbers.

Problem 1.3	 	Define a function to calculate the distance between  two points in the plane.

Problem 1.4		Define a function to compute the sum of the first n integers. Hint: Mathematica has a function called Sum.

Problem 1.5 	Plot the natural log of the quantity sin(x)+2 from x=0 to x=4π. Label your axes.

Problem 1.6		Plot tan(x) from 0 to 2π. Play with the option PlotRange.

Problem 1.7		Plot sin(x^2+y^2) with x and y running from -π to π. Play with the option  PlotPoints.

Problem 1.8		Make a density plot of sin(x^2+y^2) with x 	and y running from -π to π. Play with the option PlotPoints.

Problem 1.9 	Make a contour plot of sin(x^2+y^2) with x and y running from -π to π. Play with the option PlotPoints.

Problem 1.10 	Define a function that plots a function from 0 to 1. That is write a function called, myplot, which produces (for example)

myplot[Sin[6x+3]+x,x]

Integration

Analytic Integration

Mathematica can do both definite and indefinite integrals. The command is Integrate. The syntax for an  indefinite integral is Integrate[f[x],x] and for a  definite integral is Integrate[f[x],{x,xmin,xmax}]. For example:

Integrate[Tan[x],x]

Clear[f]
f[x_]=Log[x]/x

Integrate[f[x],x]

If Mathematica cannot do the integral it will echo it back.

Integrate[x^x,x]

Here is an example of definite integration.

Integrate[Sin[x],{x,0,Pi/2}]

AreaUnderCurve=Integrate[f[x],{x,1,2}]

N[AreaUnderCurve]

It is always a good idea to check an integral (especially a definite integral) whenever possible. This can be difficult, but the case at hand is easy to check (at least for reasonableness) just by plotting it.

Plot[f[x],{x,1,2}]

Integrals that depend on a parameter may return a conditional expression as an answer.

Integrate[Exp[-a x],{x,0,Infinity}]

This can be prevented by using the option Assumptions.

Integrate[Exp[-a x],{x,0,Infinity},Assumptions->{a>0}]

Multiple integrals can be done simply by nesting Integrate. For example

Integrate[(r^2) Sin[θ],{r,0,R},{θ,0,π}]

Using “real” mathematical notation

You may wish to use “real” mathematical notation to make your input and output easier to read. There are three different styles that you can use for input; InputForm, StandardForm and TraditionalForm. InputForm is the form that we have been using, StandardForm looks more like traditional mathematical notation and TraditionalForm looks the most like traditional mathematical notation. For output your choices are OutputForm, StandardForm and TraditionalForm. We have been using TraditionalForm for the output type. The Cell menu lets you chose between the various forms. The Convert To sub-menu allows you to convert one form to another. The Default Input FormatType and the Default Output FormatType allow you to choose the default form.

The easiest way to get (for example) an integral sign for input is to go to the Palettes menu from this menu chose Basic Math Assistant.

The next two cells are in StandardForm. Notice that I have used a proper superscript. That is, (x^2)not x^2. This is possible only in StandardForm and TraditionalForm input. Note that (on the screen) the cell brackets on the first right of the notebook have a “triangle like” mark.

∫(x^2)\[DifferentialD]x

The next two cells are in TraditionalForm. Notice the dashing on the cells brackets.

∫(x^2)\[DifferentialD]x

I usually use StandardForm for input and TraditionalForm for output. It is tempting to use TraditionalForm for input, but this is dangerous. TraditionalForm can be ambiguous. For instance is f(x), f times x or the function f(x) ?

Difficulties with Integrate

Integration can give computer algebra systems trouble and all computer algebra systems occasionally return wrong or incomplete answers. In addition Mathematica assumes that the parameters in your integral have generic values. For example

Integrate[1/x^n,x]

returns a result that is incorrect for n=1.

The following integral

result=Integrate[1/(2+Cos[x]),x]

is correct only for -π≤x≤π as can be seen by plotting  “result” along with the correct answer. The option PlotStyle (along with the graphics directives RGBColor and GrayLevel) have been used to color the incorrect curve red and the correct curve black.

correctResult=-2 (3)^(1/2) ArcTan[Sin[x]/
(Cos[x]+1)]/3+2 (3)^(1/2) ArcTan[
(3)^(1/2) Sin[x]/(3(Cos[x]+1))]/3+(3)^(1/2) x/3

Plot[{result,correctResult},{x,-2π,2π},PlotStyle->{RGBColor[1,0,0],GrayLevel[0]}]

In general, branch cuts can give Integrate trouble. When doing definite integrals you need to be careful with the integration of  functions which are continuous, but not smooth.

Numerical Integration

To do a definite integral numerically, use NIntegrate instead of Integrate. You should always plot the integrand before you do a numerical integral to see if it has singularities or if it is highly oscillatory . For example, the following integral is easy to do numerically,

Clear[f]
f[x_]=Log[x]/(1+x^2)

Plot[f[x],{x,1,2}]

NIntegrate[f[x],{x,1,2}]

but the next example is not. In version eight NIntegrate still handles the integral quickly  and accurately but this is not the case in earlier versions.

NIntegrate[Sin[1/(x^(7/2))],{x,0,1}]

Plotting shows you why.

Plot[Sin[1/x^(7/2)],{x,0,1}]

Documentation describing how NIntegrate works can be found here.

Simple Algebra

One of Mathematica's must useful features is its ability to do algebra. Below are a few of the most useful commands for doing algebra.

Factor

Factor does just what the name implies: it factors polynomials.

a=Factor[12 + 3x - 16x^2 - 16x^3 - 11x^4 + 2x^5 + 5x^6 + x^7]

Expand

Expand does just the reverse: it expands expressions.

Expand[a]

Expand[(1+x/y)^7]

Expand[(Sin[x]+Cos[y])^3]

Collect

Collect collects up terms. Its use is best understood through some examples.

b=Expand[(1+x+y^2)^3]

Collect[b,x]

Collect[b,y]

c=Expand[(x+y)^2 (Sin[x]+Cos[x])+(1-Cos[x])^3]

Collect[c,{Sin[x],Cos[x]}]

Simplify

One of the most useful of all Mathematica commands is Simplify. Simplify tries to write an expression in the simplest possible form. For example:

d=Sin[x]^2+Cos[x]^2

Simplify[d]

A more powerful form of Simplify is FullSimplify. FullSimplify is slower than Simplify, but works on special functions and applies more rules than Simplify. For instance the expression

a=((Sqrt[((3-x)/(2+x))])/(Sqrt[3-x]))

can obviously be simplified, but Simplify fails to do so.

Simplify[((Sqrt[((3-x)/(2+x))])/(Sqrt[3-x]))]

FullSimplify does the job.

FullSimplify[((Sqrt[((3-x)/(2+x))])/(Sqrt[3-x]))]

Simplify can also be used with Assumptions.

Simplify[Sin[n π],n∈Integers]

Replacement Rules

A powerful feature of Mathematica is the ability to replace a term or pattern in an expression by another expression. This is done by the command   ReplaceAll which has the short hand notation /. The syntax is expression/.old part->new part. Suppose I have an expression (say the fifth Legendre polynomial involving cos(θ) and I want to change variables to x=cos(θ).

a=LegendreP[5,Cos[θ]]

b=a/.Cos[θ]->x

The replacement is local. That is  a is unchanged.

a

b

Series Expansions

To expand a function in a Taylor series use the command Series. The syntax is Series[f[x],{x,point around which you want to expand,order}]. For example:

approximation=Series[Log[1+Sin[x]^2],{x,0,8}]

Notice that this expression has been tagged with an order of symbol. To plot this expression we must remove the order of symbol. This is done by  using Normal.

approximation=Normal[approximation]

Now we can compare the function and its series expansion. The command PlotStyle tells Mathematica  to use the first style in the list for the first function, the second style for the second function and so on. RGBColor stands for red, blue, green color. The arguments to RGBColor must be between 0 and 1.

Plot[{Log[1+Sin[x]^2],approximation},{x,0,1},
	PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}]

Exercise Set Two

Problem 2.1		Integrate sin(x)/(sin(x)+cos(x)) with respect to x.

Problem 2.2		Integrate ln(x)/(x^2) from 1 to ∞ with respect to x.

Problem 2.3 		Integrate sin((x^3)) from x=0 to x=1.

Problem 2.4		Plot (∫_0^1)sin((Sqrt[(x^(2 ))+δ]))dx with δ running from 0 to π.

Problem 2.5		Write (((x+1))^5)+((((x^2)-y))^4) in powers of x. Rewrite it in powers of y.

Problem 2.6		Write a replacement rule or rules that will rewrite Log[a] +Log[b] as Log[ab] and Log[a]-Log[b] as Log[a/b] for any a and b.

Problem 2.7		Expand Sin[(x^2)+1] about x=1. Do this for several 	orders and plot these approximations on the same graph as the function.

Differentiation

To differentiate a function you use D. If you just want  the first derivative, you can use D[f[x],x]. Thus

D[x^2,x]

D[x^3,{x,2}]

As a more complicated example let's check an integral we did earlier.

Clear[f]
f[x_]=Log[x]/x

answer=Integrate[f[x],x]

D[answer,x]

Limits

To take a limit in Mathematica you use the command Limit. The syntax for Limit is Limit[f[x],x->x0]. For example:

Limit[Sin[x]/x,x->0]

Note that Limit takes a one-sided limit! The direction is controlled by the option Direction which takes the values -1 (limit from above) or 1( limit from below). This means that Limit will return an incorrect result if the single sided limits exist but the limit does not. For example

f[x_]=x/Abs[x]

Plot[f[x],{x,-1,1}]

Limit[f[x],x->0]

Limit[f[x],x->0,Direction->-1]

Limit[f[x],x->0,Direction->1]

As a more interesting example, let us build our own derivative operator. We begin by defining the derivative of a function f[x]. We need wild card variables for both f and x in order to match whatever is in the slots.

Clear[f,h,x,d]
d[f_[x_],x_]:=Limit[(f[x+h]-f[x])/h,h->0]

Now, we need to test our operator.

d[Sin[x],x]

d[Log[x],x]

d[Tan[y],y]

As it is written, our derivative will not work on an expression like (x^2) since the operator  d expects something of the form f[x].

d[x^2,x]

Nor will it work on (e^x). To fix this need a better definition for d. The following definition provides a more general match.

Clear[d]
d[f_,x_]:=Limit[((f/.x->(x+h))-f)/h,h->0]

Now our derivative will work on (e^x) and (x^2).

d[x^2,x]

d[x^2+x,x]

d[x Exp[x^2],x]

d[Sin[x] Exp[x^2],x]

To compute a residue in Mathematica you use the command Residue. The syntax for Residue is Residue[f[z],{z,z0}] which computes the residue of f[z] at (Subscript[z, 0]). For example:

Clear[a,z,b]
Residue[z/(z-a)(z-b),{z,a}]

Solving Equations

Analytic Equation Solving

The function Solve will solve a polynomial equation of fourth order or lower, a system of linear equations, and some other equations. The syntax is Solve[equation,{variable or variables to solve for}]. You must give Solve an equation.  Mathematica makes a distinction between an assignment statement and an equation. A double == is used for an equation. Thus x=1 is an assignment of the value 1 to the variable x but x==1 is an equation whose solution is x=1. For example:

Clear[a,b,c]
solutions=Solve[a x^2+b x+c==0,x]

Notice that I get both solutions of the quadratic equation and that they are given as replacement rules. To use these solutions I have be able to pick out the one I want. To get a part of an expression in Mathematica you use expression[[part number]]. Thus solutions[[1]] gives the first solution and solutions[[2]] gives the second solution.

solutions[[1]]

solutions[[2]]

Now we need to go from a replacement rule to an assignment statement. I will call the two solutions x1 and x2.

x1=x/.solutions[[1]]

x2=x/.solutions[[2]]

To solve a system of equations you give Solve a list of equations. For example:

a=Solve[{x+y==2,x y==3},{x,y}]

x1=x/.a[[1]]

y1=y/.a[[1]]

For higher order equations Solve will generally return a Root object. For example

solutions=Solve[(x^7)+(x^3)+7==0,x]

Root object may seem strange but they just represent the roots of equations. For example Root[((#1)^2)-2,1] is the first root of the equation (x^2)-2=0 and is thus -(Sqrt[2]).

Root[((#1)^2)-2,1]

Root[((#1)^2)-2,2] is the second root.

Root[((#1)^2)-2,2]

Root objects can numericized efficiently and accurately. The roots of our seventh order polynomial are (to 30 digits)

N[solutions,30]

In version eight (but not earlier versions) Mathematica supports transcendental root object and solving equations or domains. Suppose we want to find all the roots of (Subscript[J, 1])(x)=0 on the interval 10<x<40. All we need to do is (note that && is the logical “And “operator)

besselRoots=Solve[BesselJ[1,x]==0&&10<x<40]

besselRoots//N

Numerical Equation Solving

The function FindRoot solves equations numerically. The syntax is FindRoot[equation,{x,x0}] where (Subscript[x, 0]) is the starting point for the solution search. FindRoot will find only one root at a time and may find a root other than the one you want. Here we find the first three zeros of (x^2) Sin[x] -x +1 =0.

Plot[(x^2) Sin[x]-x+1,{x,0,10}]

FindRoot[x^2 Sin[x]-x+1,{x,2.5}]

FindRoot[x^2 Sin[x]-x+1,{x,6}]

FindRoot[x^2 Sin[x]-x+1,{x,10}]

Vector Analysis

First we define a couple of vectors and a field.

Clear[a,b,c,e,f,g]
A={a,b,c};
B={e,f,g};
F[x_,y_,z_]={fx[x,y,z],fy[x,y,z],fz[x,y,z]};

We can take the cross product with the function Cross and the dot product with function Dot.

Cross[A,B]

Dot[A,B]

The functions Div, Grad Curl the Laplacian do exactly what their names imply. If we want to take the divergence of an vector field in Cartesian coordinates the syntax is Div[function,{variables}]. For example

Div[F[x,y,z],{x,y,z}]

and

Div[Curl[F[x,y,z],{x,y,z}],{x,y,z}]

If we want the Laplacian in Spherical coordinates we must specify both the variables and the coordinate system.

Laplacian[f[r,θ,ϕ],{r,θ,ϕ},"Spherical"]//Expand

If we need the Laplacian in elliptical coordinates we can use

Laplacian[f[μ,ν],{μ,ν},"Elliptic"]//Expand

We can discovver the supported coordinate systems with

CoordinateChartData[{ All,"Euclidean", 2}]

and

CoordinateChartData[{ All,"Euclidean", 3}]

Coordinate names can be looked up with

CoordinateChartData[{"Elliptic",2}, "StandardCoordinateNames"]

and

CoordinateChartData[{"Elliptic",2}, "AlternateCoordinateNames"]

Using Packages

Mathematica has a number of packages which contain code for various  specialized applications. The package VectorAnalysis for example, contains functions for doing vector calculus. Much of the functionality in the packages will end up in the kernel in future versions of Mathematica.

To load a package you type <<”packagename`”. In this case the directory is Calculus and the package name is VectorAnalysis.  Needs["packagename`"]. The advantage of Needs as opposed to << is that Needs will load a package only if it has not already been loaded. Thus, Needs will not load a package twice. This is important because many Mathematica packages are not designed to be read in twice. If there is a package that you use frequently, you can put the command to load it in your init.m file and it will be loaded each time you start Mathematica.

Some Mathematica  packages need to be loaded as the first line of code run in a Mathematica session. If you fail to do this, you may get error messages saying that variables are shadowed in other contexts. This simply means that you have  already used a name that is exported by the package. Your definition of the symbol will shadow the definition of the symbol given in the package and the package will not work correctly. If you have only a few shadowed symbols you can use Remove to remove the offending symbols. If you have a great many shadowed symbols it may be easier to quit and restart the kernel.  Remove[symbol] removes "symbol" from the Symbol table. If we mention a symbol Mathematica will put it in the Symbol table. Lets create a symbol called xx:

xx

Typing ?xx shows that Mathematica knows about the existence of this symbol and that xx is defined in the global context. Contexts are similar to file directories on a UNIX system.

?xx

We can remove this symbol by using Remove.

Remove[xx]

?xx

If a package loads correctly, you will not get any output from loading it.

Needs["VariationalMethods`"]

We now have access to the functions in VariationalMethods. For example the Lagrangian for a pendulum is

ℒ=(1/2)m (r^2)((D[θ[t],t])^2)+m g r Cos[θ[t]]

equationOfMotion=EulerEquations[ℒ,θ[t],t]

Which we can now use DSolve to solve.

DSolve[equationOfMotion,θ[t],t]

We can get a list of all the functions in the package by asking for all the symbols that are in the package context.

?VariationalMethods`*

Animation and Graphics Primitives and Manipulate

Animation is easy in Mathematica. You make a table of graphs and then use ListAnimate. Here we turn the sine into a cosine and back into a sine again.

Animation and Graphics Primitives

ListAnimate[Table[Plot[{Sin[x+delta],Cos[x]},{x,0,2Pi},
PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}],
		{delta,0,2Pi,Pi/8}]]

You can control the speed of the animation by using the keys on the number pad or by using the VCR controls on the graphic. This works for any set of  Mathematica graphics.

Mathematica also has graphics primitives with which you can draw. They include Point, Line, Circle, Rectangle, Cuboid, Polygon and Disk among others. To use them you need to make the graphics primitive a graphics object.  You need to use  Graphics on two-dimensional objects and Graphics3D on three dimensional objects in order to make them graphics objects.

object={Line[{{0,0,0},{1,0,0}}],Line[{{0,0,0},{0,1,1}}]}

lines=Graphics3D[object]

point=Graphics3D[{Red,Point[{1/2,1/2,1/2}]}]

You can combine graphics objects using Show.

Show[point,lines,Boxed->True]

You can use RGBColor,Thickness, GrayLevel, Hue and PointSize to control the color and size of the graphics primitives.
 
 As a more realistic example, we will animate a box undergoing simple harmonic motion. We will need to make four objects. A box, a spring, a surface for the box to slide on and a wall to anchor the spring to. We will start with a box which must oscillate with time. The left edge of the box is at x=1 at t=0 and its right edge is at x=2 at t=0. It has a height of 1/2 and I have colored it red using RGBColor. Style directives such as RGBColor must come before the graphics primitive. The  style directives and the graphics primitives are then passed to Graphics or Graphics3D  (to make them graphics objects) enclosed in a list.

box[t_]:=Graphics[{RGBColor[1,0,0],
	Rectangle[{1+Sin[t],0},{2+Sin[t],1/2}]}]

Show[box[0]]

Next, we make our surface. I have put the surface a little bit below zero so that when we are finished our box will ride on it. I determined how far below zero I wanted the surface by trial and error.

surface=Graphics[{Thickness[0.05],RGBColor[0,1,0],Line[{{0,-1/20},{3,-1/20}}]}]

The option AspectRatio->Automatic has been used to ensure that the true (physical) aspect ratio is used.

Show[surface,box[1/2],AspectRatio->Automatic]

The next step is the attach our spring to the box.

spring[t_]:=Graphics[{Thickness[0.015],
	Line[{{0,1/4},{1+Sin[t],1/4}}]}]

Show[surface,box[Pi/2],spring[Pi/2],
AspectRatio->Automatic]

Lastly, we make the wall and put everything together.

wall=Graphics[{Thickness[.05],RGBColor[0,0,1],
	Line[{{0,-1/20},{0,1}}]}]

Show[surface,box[Pi/2],spring[Pi/2],wall,
AspectRatio->Automatic]

Now for the animation!

ListAnimate[Table[Show[surface,box[t],spring[t],wall,
AspectRatio->Automatic],
	{t,0,2Pi,Pi/8}]]

Manipulate

A very useful function for visualization is Manipulate. The basic syntax is Manipulate[expression,{u,umin,umax}] which evaluates “expression” with the parameter u running from umin to umax. Manipulate[expression,{u,umin,umax,step}] evaluates “expression” with the parameter u running from umin to umax in steps of “step”. A few examples should make the use of Manipulate clear. Just drag the slider bar to change the value of the parameter.

Manipulate[Expand[(((1+x))^n)],{n,1,10,1}]

Manipulate[Plot[BesselJ[μ,x],{x,-3,10},PlotLabel->StringJoin["plot of ", ToString[TraditionalForm[(Subscript[J, μ])]]]],{μ,0,5}]

A much more complicated example

The dynamic elements in Mathematica allow you to build complicated interactive graphics with out much trouble. In the example below I have plotted the field lines from three charges. You can move the gray Locator points around with your mouse. This will change the locations of the charges. Moving the slider bars will change the values of the charges. Drawing stream plots is computationally intensive so you will need to wait for a second or so for the plot to re-draw.

Manipulate[DynamicModule[{locations={{1,0},{-1,0},{0,1}}},
LocatorPane[Dynamic[locations],
Dynamic[charges=Map[Plus[{x,y},-#]&,locations];
f=Drop[Grad[Apply[Plus,{q1,q2,q3}/Simplify[Map[EuclideanDistance[{0,0},#]&,charges],x∈Reals&&y∈Reals]],{x,y,z}],-1];StreamPlot[Evaluate[f],{x,-2,2},{y,-2,2}]]]],{{q1,1},-2,2},{{q2,1},-2,2},{{q3,-1},-2,2},SaveDefinitions->True,ContinuousAction->False,SynchronousUpdating->False]

ss

Solving Differential Equations

Analytically

The function  DSolve solves differential equations analytically. There is of course, no guarantee that  DSolve will be able to solve your equation. The syntax for DSolve is DSolve[{eqn1,eqn2,....}, {y1[x1,x2,...],y2[x1,x2,...],...}, {x1,x2,...}]. If you have only one equation, the form is  DSolve[equation,y[x],x] or DSolve[{equation,initialconditions},y[x],x], depending on whether you give the initial conditions.

Suppose we have a simple harmonic oscillator. The equation we want to solve is:

SHO=D[x[t],{t,2}]+ω x[t]==0

answer=DSolve[SHO,x[t],t]

The constants c1 and c2 are determined by the initial conditions. Like all of Mathematica equation solvers,  DSolve returns a replacement rule. Furthermore it returns the replacement rule as a nested list. We can get just a single list by asking for the  First part of answer.

First[answer]

x[t_]=x[t]/.First[answer]

The velocity is

Clear[v]
v[t_]=D[x[t],t]

Notice that we used wild-card variables in the definition of x[t] and v[t]. For our initial conditions, suppose we start the oscillator at x[0]=1 and v[0]=0.

ICs=Solve[{x[0]==1,v[0]==0},{C[1],C[2]}]

x[t_]=x[t]/.ICs

v[t_]=v[t]/.ICs

A simpler way of dealing with the ICs is to include them in the call to DSolve.

Clear[x,v]
answer=DSolve[{SHO,x[0]==1,x'[0]==0},x[t],t]

We could also treat this as a system of equations.

Clear[v,x]
answer=DSolve[{v'[t]+ω x[t]==0,v[t]==x'[t],x[0]==1,x'[0]==0},{v[t],x[t]},t]

Series Solutions

Consider once again the simple harmonic oscillator (((d^2)y)/(dx^2)) + (ω^2) y =0 . Let us find a series solution to this equation.

Clear[eqn,y,x,ω,k,F,a]
	eqn=D[y[x],{x,2}]+(ω^2) y[x]

We seek a solution of the form y(x) = (x^k) (∑_(n=0)^∞)(Subscript[a, n]) (x^n). Expanding y(x), to n=5 gives

y[x_]=(x^k) (Sum[a[n] (x^n),{n,0,5}]+((O[x])^6))

and collecting terms in x gives

Collect[eqn,x]

shows us that the indicial equation is k(k-1)=0 and thus we must have either k=1 or k=0. We turn first to the k=1 case.

k=1

The function LogicalExpand can be used to equate coefficients in x.

equations=LogicalExpand[eqn==0]

The ∧ symbol is the logical and ( you could also use &&).

These equations can be solved using Solve. We want to solve for all the coefficients except a[0] and we can use Table to make a list of the variables we want to solve for.

variables=Table[a[n],{n,1,5}]

results=Solve[equations,variables]

We can use ReplaceAll (also known as /.) to go from a replacement rule to a assignment statement.

y[x]/.results

Removing a set of curly brackets by using First and getting rid of the order of symbol by using Normal gives

firstSolution=Normal[First[y[x]/.results]]

If we set ω and a(0) to 1 we can compare our solution to the exact solution.

Plot[{firstSolution/.{a[0]->1,ω->1},Sin[x]},{x,0,2Pi},PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}]

If we put all the code together into a single cell and use m to denote the order of the series expansion we can change the order by changing a single number. I have used semi-colons to suppress all but the last line of output and set m to 13.

Clear[eqn,y,x,ω,k,F,m,firstSolution]
eqn=D[y[x],{x,2}]+(ω^2) y[x];
m=13;
y[x_]=(x^k) Sum[a[n] (x^n),{n,0,m}]+((O[x])^(m+1));
k=1;
equations=LogicalExpand[eqn==0];
variables=Table[a[n],{n,1,m}];
results=Solve[equations,variables];
firstSolution=Normal[First[y[x]/.results]];
Plot[{firstSolution/.{a[0]->1,ω->1},Sin[x]},{x,0,2Pi},PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}]

The k=0 case works in almost exactly the same way. Since we have a two term recursion relationship that relates (Subscript[a, n+2 ])to (Subscript[a, n]) , the a[1] term is arbitrary (this  of course was true for the k=1 case as well) and we will set it to zero. Making this change, setting  k=1 and m=12 gives

Clear[eqn,y,x,ω,k,F,m,secondSolution,a]
eqn=D[y[x],{x,2}]+(ω^2) y[x];
m=12;
y[x_]=(x^k) Sum[a[n] (x^n),{n,0,m}]+((O[x])^(m+1));
k=0;
equations=LogicalExpand[(eqn/.a[1]->0)==0];
variables=Table[a[n],{n,2,m}];
results=Solve[equations,variables];
secondSolution=Normal[First[y[x]/.results]/.a[1]->0];
Plot[{secondSolution/.{a[0]->1,ω->1},Cos[x]},{x,0,2Pi},PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}]

Numerically

Mathematica can solve certain classes of differential equations symbolically but many, in fact most, differential equations have no analytic solution and must be attacked either numerically or by using some sort of approximation scheme. We will consider Duffing's equation. The canonical form for Duffing's equation is :

(Overscript[x, ¨])+2γ (OverDot[x])+((Subscript[ω, 0])^2) x+β (x^3) = f cos(ω t)

Physically this corresponds to a damped non-linear oscillator with a driving force. The spring force is

F(x) =-m(((Subscript[ω, 0])^2)x + β (x^3))

If (Subscript[ω, 0]) is real and β is positive, Duffing's equation can be put into a simpler form by rescaling the coordinates and the time. In this form Duffing's equation is:

(Overscript[x, ¨])+2γ (OverDot[x])+ x+(x^3) = f cos(ω t)

We can solve this equation by using NDSolve. The syntax is

??NDSolve

Mathematica's description of  NDSolve is not complete. You will notice that there is no mention of initial conditions, yet you cannot find a numerical solution to a differential equation without giving the initial conditions. The correct form is NDSolve[{eqns,initialdata},y, {x,xmin,xmax}].

First we will define an expression called duffing. When this expression is set equal to zero we get Duffing's equation.

duffing=D[x[t],{t,2}]+2γ D[x[t],t]+x[t]+((x[t])^3)-f Cos[ω t]

We will use NDSolve to solve this equation. The form must be eqn==0. Remember, the == tells Mathematica that this is an equation. The command duffing==0/.{γ->1/10,f->1,ω->4} tells Mathematica to replace the variables γ, f and ω by 1/10,1, and 4 respectively. The replacement occurs only in this one expression, duffing is not globally changed. Ask for "duffing" and see for yourself. The initial conditions must be given in the form of equations using ==. Notice that   Mathematica understands the prime notation for derivatives. In this example we start the particle at the origin with zero initial velocity and integrate the equation for the first 8 seconds.

ans=NDSolve[{duffing==0/.{γ->1/10,f->1,ω->4},x'[0]==0,x[0]==0},x,{t,0,8}]

The answer that we get is a replacement rule in which x has been replaced by an interpolating function. The <> in the answer tells use that there is hidden information, which is the actual interpolation function and the (0.,8.) piece tells us that the interpolation function starts a 0. and ends at 8. Let's plot the answer. To do this we replace x[t] (remember that x[t] and x are not the same) by the answer, which we have called "ans". Before plotting we  evaluate the result. The use of Evaluate speeds up the plotting.

Plot[Evaluate[x[t]/.ans],{t,0,8}]

In the example below, the value of f, the driving force, has been changed and we integrate the equation for the first 70 seconds. We then plot it (but only for the first 12 seconds) just as before.

Clear[x]
ans1=NDSolve[{duffing==0/.{γ->1/10,f->3,ω->4},x'[0]==0,x[0]==0},x,{t,0,70}]

Plot[Evaluate[x[t]/.ans1],{t,0,12},AxesLabel->{"t","x[t]"}]

We should make a phase space plot so that we  get a better feeling for the nature of the motion. To make a phase space plot, we plot v versus x (really p versus but the difference is just a scale factor). To do this we replace x[t] by the interpolating function and take the derivative. This gives us v, which is plotted  below.

v=Evaluate[D[x[t]/.ans1,t]]

Plot[v,{t,0,12},AxesLabel->{"t","x[t]"}]

We will also need x.

x=Evaluate[x[t]/.ans1]

To make our phase space plot we want to make a parametric plot of v(t) versus x(t). The Mathematica command for this is  ParametricPlot.  The syntax is ParametricPlot[{fx,fy},{t,tmin,tmax}]. There is one slight complication. If we ask for v at some time, say one second,we get.

v/.t->1

The answer has braces around it.  ParametricPlot wants a number with no braces! We can get just the number by asking for the first part of v.

v[[1]]/.t->1

ParametricPlot[{x[[1]],v[[1]]},{t,0,70},
AxesLabel->{"x[t]","v[t]"}]

This phase space plot looks non-periodic at first but you can see that the later orbits settle down and the phase space plot looks a lot like a linear oscillator. For larger f the phase space plot is more complicated as we see in the next examples.

In the example below  all the commands  are combined in a single cell. We must first clear x, x[t] and x'[t], otherwise NDSolve thinks that they are interpolating functions. Clearing x will also clear x[t] and x'[t]. In order to wash out the initial conditions, and look only at the asymptotic behavior of the orbits, I have plotted the motion for the range 40 to 70 seconds.  Bundling the code together into one gives

Clear[x]
ans2=NDSolve[{duffing==0/.{γ->1/10,f->20,ω->1.45},x'[0]==0,x[0]==0},x,{t,0,70}];
v=Evaluate[D[x[t]/.ans2,t]];
x=Evaluate[x[t]/.ans2];
ParametricPlot[{x[[1]],v[[1]]},{t,40,70},PlotPoints->200]

Here we start to see period doubling. Period doubling is said to occur when it takes two cycles for the orbit to close on itself. In the plot above you see two dark lines which is the signal for period doubling.

The phase space plots for large f are typically very complex as is shown in the example below. I have used the option MaxSteps to

Clear[x]
ans2=NDSolve[{duffing==0/.{γ->1/10,f->100,ω->1.40},x'[0]==0,x[0]==0},x,{t,0,120},MaxSteps->20000];
v=Evaluate[D[x[t]/.ans2,t]];
x=Evaluate[x[t]/.ans2];
ParametricPlot[{x[[1]],v[[1]]},{t,70,120},PlotPoints->200]

Solving a System of Equations

You can use NDSolve to solve a system of differential equations. For a system of equations the syntax is:
 NDSolve[{eqn1==0,eqn2==0,...,eqn3==0,intialconditions},
 {y1,y2,...yn},{x,xmin,xmax}].

Exercise Set Three

Problem 3.1 	Check the result for the integral of sin(x)/(sin(x)+cos(x)).

Problem 3.2 	Compute the limit of (x^2) csc(x) as x approaches 0. Compute the limit of (x-1)/(Sqrt[(((x-1))^2)]) as x approaches 1.

Problem 3.3	Find  all the roots of sin((x^2))-cos(x) on the interval 0 to π.

Problem 3.4	 Solve  the differential equation: 	
		f''(x)-4f'(x)+4f(x)=0 with f(0)=1 and f'(0)=1. Plot the result.

Problem 3.5	Solve  the differential equation: 	
		f''(t)-4f'(((t))^2)+4f(t)=0 with f(0)=1 and 	f'(0)=1 for t between 0 and 20. Plot the result.

Problem 3.6	Modify the code for the simple harmonic oscillator so that the motion of the block is governed by the solution to problem 3.5.