Slightly Advanced Introduction to Mathematica.

Below is a slightly advanced introduction to some concepts in Mathematica.  You may be familiar with most of them.

On the right hand side there are a number of brackets,  Click on the one that is beside the expression "2+2".  Then press the return key and the shift key simultaneously.  This will start  up the Mathematica "kernel" or central processing program.  It  may take a few seconds to respond.

2+2

Mathematica can do more than evaluate 2+2.  Click on the bracket to the right of each of the following expressions and  again press "shift-return".

4 (3^6)-28 (7^2)

((234 128-48 632)/(54^2))

((Sin[(π/6)])^4)+((Cos[(π/4)])^2)

100!

(Sqrt[15241578750190521])

In Mathematica adjacent terms are multiplied (you don't have to explicitly put in a multiplication symbol).  To raise something to a power, you can either use the caret "^" symbol, or use the BasicInput palette and typeset an exponent.  (Select the "Palette" submenu from the "File" menu).   Thus:

                            x^2 = (x^2)

   Also note that we have functions like Sqrt[...] for taking the square root, and sines and cosines are Sin[] and Cos[].  The capital letters are important, as are the square brackets.  All predefined functions in Mathematica are capitalized, and take square brackets.  The curly brackets {} are used to separate  parameters, or to mark elements as belonging to a list.  The  program will whine petulantly and refuse to work if you use  the wrong brackets in the wrong place, or forget a comma when  it wants one there.

The most important point in this whole presentation:

One extremely important point: highlight the word "Sin" above in the input box.  Now pull down the help menu and select "Find Selected Function".  This will tell you all about the sine function in Mathematica.  The help browser in Mathematica is extensive and extremely helpful.  You can look up functions alphabetically in its index or browse them by subject. You can take a tour of Mathematica and see what it can do.  You can see thousands of examples.  Learn to use the browser.

Mathematica is a Symbolic manipulator:

Mathematica is a symbolic manipulator.  It does not always evaluate things in a way you might expect!  Below are two  examples.  In each case we as Mathematica to evaluate an  expression.  If we want a decimal or NUMERICAL answer, we have  to use the  function N[...], which calculates its argumument in a numerical fashion.

(14/6)

N[(14/6)]

(Sqrt[68])

N[(Sqrt[68])]

(Sqrt[68.])

Mathematica is more than just a fancy calculator.  It can perform algebraic operations as well.  The form or "syntax"   is:

Solve[{eq1==rhs1,eq2==rhs2,...},{var1,var2..}]

Note the funny double equal signs (==).  This is used because we are specifying a symbolic relationship, we're not just renaming the eq1 by rhs1.   There are a few examples below.  Again, click on the bracket  on the RHS and hit shift-enter.
  
  A quadratic equation:

Solve[{(x^2)+x-12==0},{x}]

A cubic equation:

Solve[{(x^3)-3 (x^2)+x+1==0},{x}]

Even a quartic equation!

Solve[{(x^4)+(x^3)-4 (x^2)-x+1==0},{x}]

You can also solve systems of equations this way:

Solve[{x+y == 12, 4x-2y== 24},{x, y}]

However this is not the recommended way for large systems of equations - use

Integrals:

Mathematica can do integrals symbolically.   Here's an integral you probably don't want to do by hand:

Integrate[Sin[x]/x,{x,0,Infinity}]

Mathematica can do this integral numerically.

NIntegrate[Sin[x]/x,{x,0,Infinity}]

The point is that Mathematica can also do integrals that have no hope of being done analytically:

NIntegrate[((BesselJ[0,x])^2) Sqrt[(1/(1+(x^2)))] ,{x,0,1}]

Vectors and Matrices:

The fundamental item in Mathematica is the "list."   This is much like a vector. It consists of a set of elements  enclosed in curly brackets:  "{a1, a2, a3...}".    We can construct a list either by hand, or by having Mathematica perform some action repeatedly using the Table command.  The latter is much like a "do-loop" in Fortran or a "for-next" loop in Basic.

vector1={1,2,3,4}
vector2=Table[n!,{n,1,20}]

When working with vectors and matrices, it is important to know when to tell Mathematica to shut up.  If you put a semicolon on the end of a line, then the calculation is performed, but the output is suppressed.

vector1={1,2,3,4};
vector2=Table[n!,{n,1,20}];

If you want a pretty version of the vector, then use MatrixForm.

MatrixForm[vector2]

To refer to an element in the list, we just double square brackets.  For example, the 10th element in vector2 is vector2[[10]]:

vector2[[10]]

A matrix is just a list of lists.  Again, we can write it out by hand or have Mathematica create it using a Table command.

matrix1={{1,2,3,4},{5,0,7,8},{5,4,0,2},{7,6,5,4}};
matrix2= Table[  Sin[i*j],{i,1,4},{j,1,4}];
MatrixForm[matrix1]
MatrixForm[matrix2]

Again, we can address an individual element using the double square brackets, but now we have to have two indices:

matrix1[[2,3]]

Mathematica can handle lots of matrix manipulations.  Here’s the inverse of a matrix:

Inverse[matrix1]

Mathematica will find exact solutions, often when we don't really want them.  Please note the difference between the two commands below and in their output.  Students often write slow Mathematica code because they forget to tell Mathematica to work numerically and not symbolically.

Inverse[matrix2]
Inverse[N[matrix2]]

Mathematica can solve linear sets of equations, and quickly find eigenvalues and eigenvectors.

LinearSolve[matrix1,vector1]

How fast is Mathematica?   Let’s solve full matrices of a size (2^n)*(2^n) in size:

nval=4;
bigMatrix=Table[RandomReal[],{i,1,(2^nval)},{j,1,(2^nval)}];
bigVector=Table[RandomReal[],{i,1,(2^nval)}];
bigAnswer=LinearSolve[bigMatrix,bigVector];

Let’s see how this scales:

timeRuns=Table[bigMatrix=Table[RandomReal[],{i,1,(2^n)},{j,1,(2^n)}];
bigVector=Table[RandomReal[],{i,1,(2^n)}];
{n,Part[Timing[bigAnswer=LinearSolve[bigMatrix,bigVector];],1]},{n,1,12}]

ListLogPlot[timeRuns,PlotRange-> All]

matrix3=N[matrix1];
Eigenvalues[matrix3]
Eigenvectors[matrix3]

Are the eigenvectors normalized?  Let's check:

eigvec=Eigenvectors[matrix3];
Table[eigvec[[i]].Conjugate[eigvec[[i]] ],{i,1,4}]

Pretty darn close!  Let's chop off anything less that (10^(-10)).

Table[Chop[eigvec[[i]].Conjugate[eigvec[[i]] ]],{i,1,4}]

Mathematica cares about Capitalization

Every function or command in Mathematica begins with a Capital Letter.  Any English word in the command will have an internal Capital Letter.  The function for a contour plot is ContourPlot, not Contourplot.  Mathematica will interpret any incorrect spelling or capitalization as a definition of some new command.  For example, only one of the following is correct.  Which one?

arcsin[1]
Arcsin[1]
ArcSin[1]

Mathematica is a symbolic manipulator.

What is the value of Sin[(π/4)] ?

a=Sin[(π/4)]

Mathematica will keep this symbolic expression, unless you tell it to evaluate it as a number:

N[a]

The two expressions below are not equivalent to Mathematica: one is a collection of symbols, and one is a number.

10/7
10./7

Mathematica never forgets

Mathematica remembers everything you've executed.  Even if you've deleted the cell or forgotten about the command, it will remember it (until you either restart the kernel or tell it explicitly to forget something).   Now, let's say you want to solve the quadratic equation:

Solve[{a (x^2)+b x+c ==0},{x}]

Oops!  We don't get the quadratic formula.  We set a numerical value to "a" in the previous section.  We have to clear the value of "a"

Clear[a];
Solve[{a (x^2)+b x+c ==0},{x}]

Assignment versus equals

Look at the two expressions below.  The first assigns the value "0" to the variable "zero".  Whenever you type the variable name, Mathematica just substitutes in the value you assigned.  The second expression is an equation.  It is a logical relation.  It's not a rule for substitution, but a relationship that you can use in solving systems of equation.

zero=0
(y1^2) + 3 == 0

Plotting against you

Mathematica can only plot numbers.  If any part of your expression is not defined, it won't be able to plot it.

y=A Sin[k x-w t];
A=2;
k=π;
w=2;
Plot[y,{x,0,2}]

Set the value of t above and then try it again.

Variable blindness

If you don't leave spaces between your variables, Mathematica just assumes that this is some new variable.  Note that "mg" is not defined, but "m g" is defined.

m=2;
g = 9.8;
mg
m g

Combining Plots

The simplest way to combine 2D plots is to give them as a list to the plot command.

Plot[{x,(x^2),(x^3),(x^4)},{x,0,1}]

If you name individual plots, you can also combine them using Show.

plt1=Plot[x,{x,0,1}];
plt2=Plot[(x^2),{x,0,1}];
Show[plt1,plt2]

plt3=Plot3D[Sin[x y],{x,-π,π},{y,-π,π},PlotPoints->40]
plt4=Plot3D[Cos[x y],{x,-π,π},{y,-π,π},PlotPoints->40]
Show[plt3,plt4]

Differential Equations

Simple ODE’s

eqn=D[x[t],t,t] + x[t] == 0
DSolve[eqn,x[t],t]

eqns={D[x[t],t,t] + x[t] == 0,x[0]==0, x[1]==1};
DSolve[eqns,x[t],t]

Harder ODE’s

eqns={D[x[t],t,t]+ Cos[t] x[t]== 0, x[0]==0, x'[0]==1}
soln=DSolve[eqns,x[t],{t,0,16}]

Plot[Evaluate[x[t]/.soln],{t,0,16},PlotRange->All]

Non-linear DE’s

eqns={D[x[t],t,t]+ Cos[x[t]+t] +  t x[t] == 0, x[0]==0, x'[0]==1}
soln=NDSolve[eqns,x[t],{t,0,100}]
Plot[Evaluate[x[t]/.soln],{t,0,20},PlotRange->All]

Solving a discretized problem:

If we want to solve the problem discussed in class, the electron in a harmonic oscillator, we need to construct the matrix, find the eigenvalues, and look at the results.

npts=20;
dx=(1./(npts+1)); 
hmat=Table[0,{i,1,npts},{j,1,npts}];
Do[hmat[[i,i]] = kappa(((i dx -.5))^2)  +(2/(dx^2)),{i,1,npts}]
Do[hmat[[i,i+1]] =-(1/(dx^2)),{i,1,npts-1}]
Do[hmat[[i+1,i]] =-(1/(dx^2)),{i,1,npts-1}]

MatrixForm[hmat]

hmatVal=hmat/. {kappa-> 0};

Eigenvalues[hmatVal]/(π^2)

npts=200;
dx=(1./(npts+1)); 
hmat=Table[0,{i,1,npts},{j,1,npts}];
Do[hmat[[i,i]] = kappa(((i dx -.5))^2)  +(2/(dx^2)),{i,1,npts}]
Do[hmat[[i,i+1]] =-(1/(dx^2)),{i,1,npts-1}]
Do[hmat[[i+1,i]] =-(1/(dx^2)),{i,1,npts-1}]

hmatVal=hmat/. {kappa-> 0};
Eigenvalues[hmatVal]/(π^2)

evecs=Eigenvectors[hmatVal];

ListPlot[{evecs[[npts]],evecs[[npts-1]]}]

Turn on the spring:

npts=200;
dx=(1./(npts+1)); 
hmat=Table[0,{i,1,npts},{j,1,npts}];
Do[hmat[[i,i]] = kappa(((i dx -.5))^2)  +(2/(dx^2)),{i,1,npts}]
Do[hmat[[i,i+1]] =-(1/(dx^2)),{i,1,npts-1}]
Do[hmat[[i+1,i]] =-(1/(dx^2)),{i,1,npts-1}]

hmatVal=hmat/. {kappa-> 10000};
evals=Eigenvalues[hmatVal]

Let’s see if the scaling we thought was correct for the eigenvalues is actually correct.

evals/(Sqrt[10000])

Now we can plot the lowest two energy eigenvectors.   Remember that Mathematica calculates them in order of decreasing absolute value of the eigenvalues.  Therefore our last two eigenvectors are the lowest two eigenstates.

evecs=Eigenvectors[hmatVal];

ListPlot[{evecs[[npts]],evecs[[npts-1]]}]

Finally, let’s sweep the spring constant and see when we switch over from infinite square well energies to harmonic oscillator energies.  First we will calculate all of the eigenvalues.  Then we’ll pick out the ones we want.

maxn=15;
allEvals=Table[hmatVal=(hmat/. kappa-> (2^n));
Eigenvalues[hmatVal],{n,0,maxn}];

list1=Table[{(2^n),allEvals[[n+1,npts]]},{n,1,maxn}];
list2=Table[{(2^n),allEvals[[n+1,npts-1]]},{n,1,maxn}];
ListLogLogPlot[{list1,list2}]

Ok, we see that it switches behavior around κ = 100.  Can we see how it depends on κ?  First we need to take the logarithm of the x and y points.

logList1=Log[list1];
ListPlot[logList1]

Next, we need to get the last few of them, where it’s linear.

data=Take[logList1,-6]

Now we fit is to a line:

Fit[data,{1,x},x]

Solving an iterative problem:

Recall our discusion of solving the 1D problem:

	y’ = - c y  
	y(0) = 1
	
The Euler method yields:

	(Subscript[y, i+1])= (Subscript[y, i]) - c Δt (Subscript[y , i])  
Let’s test this with c = 4 and Δt = .1   

Below I iterate the problem forward in time, saving the answers in a list, euler.  I then generate plots for the list and the exact answer,  Exp[- c t], but I don’t display either plot.  After that, I use the Show command to display both plots in the same frame.

yold=1.;
c = 5.;
npts=50;
dt = 2.0/npts;
euler=Table[ynew= (1- c dt) yold ;  
yold= ynew; {i dt, yold},{i,1,npts}];
eulerPlot=ListPlot[euler,PlotRange-> All];
exactPlot=Plot[Exp[- c t],{t,0,2},PlotRange-> All];
Show[eulerPlot,exactPlot,PlotRange-> All]

Try running the routines below, changing npts to lower and lower values.

yold=1.;
c = 5.;
npts=10;
dt = 2.0/npts;
euler=Table[ynew= (1- c dt) yold ;  
yold= ynew; {i dt, yold},{i,1,npts}];
eulerPlot=ListPlot[euler,PlotRange-> All];
exactPlot=Plot[Exp[- c t],{t,0,2},PlotRange-> All];
Show[eulerPlot,exactPlot,PlotRange-> All]

Now let’s try an implicit scheme.  Note that it is always stable.  However, it is not always accurate.

yold=1.;
c = 5.;
npts=10;
dt = 2.0/npts;
euler=Table[ynew= yold/(1+ c dt) ;  
yold= ynew; {i dt, yold},{i,1,npts}];
eulerPlot=ListPlot[euler,PlotRange-> All];
exactPlot=Plot[Exp[- c t],{t,0,2},PlotRange-> All];
Show[eulerPlot,exactPlot,PlotRange-> All]

yold=1.;
c = 5.;
npts=10;
dt = 2.0/npts;
euler=Table[ynew= yold (1- (c/2) dt)/(1+ (c/2) dt) ;  
yold= ynew; {i dt, yold},{i,1,npts}];
eulerPlot=ListPlot[euler,PlotRange-> All];
exactPlot=Plot[Exp[- c t],{t,0,2},PlotRange-> All];
Show[eulerPlot,exactPlot,PlotRange-> All]

Let’s compare them all:

yold=1.;
c = 5.;
npts=20;
dt = 2.0/npts;
eulerForward=Table[ynew= (1- c dt) yold ;  
yold= ynew; {i dt, yold},{i,1,npts}];
yold=1.;
eulerBackward=Table[ynew= yold/(1+ c dt) ;  
yold= ynew; {i dt, yold},{i,1,npts}];
yold=1.;
eulerMidpoint=Table[ynew= yold (1- (c/2) dt)/(1+ (c/2) dt) ;  
yold= ynew; {i dt, yold},{i,1,npts}];
eulerPlot=ListPlot[{eulerForward, eulerBackward,eulerMidpoint},PlotRange-> All];
exactPlot=Plot[Exp[- c t],{t,0,2},PlotRange-> All];
Show[eulerPlot,exactPlot,PlotRange-> All]

Why does the midpoint approach work the best?    Well, in each time step we are approximating (e^(-c dt)).   The exact answer is:

Clear[c, dt];
Series[Exp[-c dt],{dt, 0, 6}]

Let's look at how well they do. The forward Euler approximates it by (1 - c*dt).   It has no higher order terms.  The backward approximates it by:

Series[(1/(1+c dt)),{dt,0,6}]

The above series has the (dt^2) term, but too much of it.

Series[(1- (c/2) dt)/(1+ (c/2) dt), {dt,0,6}]

The midpoint form has exactly the right amount to order (dt^2)