Intermediate Mathematica

Author: Barry McQuarrie
University of Minnesota, Morris
Jan 2014

This document provides some more advanced examples of how to use Mathematica that come up in courses beyond Calculus (in my world, that means differential equations, calculus III, and numerical methods). However, you may have some previous experience with the syntax discussed here.
 
If you have never used Mathematica before, you should begin with the Mathematica Basics file.

Partial Derivatives

Here is a good way to do partial derivatives:

Derivative[1,0][f][t,y]

Derivative[0,1][f][t,y]

f[t_,x1_,x2_,x3_]=Exp[Cos[x1*Sqrt[x2-x3]]]*BesselJ[0,t]

Derivative[0,0,0,1][f][t,x1,x2,x3]

Derivative[1,1,1,1][f][t,x,y,z]

Direction and Vector Fields

Here we define a direction field:

L[y_]:=D[y,{x,1}]==((Cos[x])/(Exp[-(y^2)]))

L[y[x]]

The direction field plots small arrows with slope give by <Δx,Δy>~<dx,dy>.

fieldeqn=VectorPlot[{Cos[x],Exp[-(y^2)]},{x,-2,2},{y,-2,2},ImageSize->{500,500}]

The size of the arrows represent the strength of the vector field. If you want all the lines to be the same length, you can use the following (which is usually a more useful plot):

fieldeqn=VectorPlot[{Cos[x],Exp[-(y^2)]},{x,-2,2},{y,-2,2},ImageSize->{500,500},VectorScale->{Tiny,Small,None}]

Differential Equations

Solving Differential Equation using DSolve

DSolve[4*y''[t]+2*y'[t]+y[t]==Sin[t],y[t],t]

Solving Initial Value Problems using DSolve

Just include the initial conditions (and add some curly brackets) if you have an initial value problem.

DSolve[{4*y''[t]-2*y'[t]+y[t]==Sin[t],y[0]==1,y'[0]==2},y[t],t]

You can then cut and past to define the solution if you want to work with it some more.

sol[t_]=(1/13) (11 (E^(t/4)) Cos[(((Sqrt[3]) t)/4)]+2 Cos[t] ((Cos[(((Sqrt[3]) t)/4)])^2)-3 ((Cos[(((Sqrt[3]) t)/4)])^2) Sin[t]+35 (Sqrt[3]) (E^(t/4)) Sin[(((Sqrt[3]) t)/4)]+2 Cos[t] ((Sin[(((Sqrt[3]) t)/4)])^2)-3 Sin[t] ((Sin[(((Sqrt[3]) t)/4)])^2))

Plot[sol[t],{t,0,20}]

Solving Differential Equation Numerically using NDSolve

Here is a nonlinear differential equations Mathematica cannot solve:

DSolve[{y''[t]+ y[t]^3==3Cos[t],y[0]==2,y'[0]==3},y[t],t]

You can get a numerical solution in an interval of your choosing (I chose 0≤t≤100). Note the format is slightly different than DSolve.

NDSolve[{y''[t]+ y[t]^3==3Cos[t],y[0]==2,y'[0]==3},y,{t,0,100}]

The following form easily gets you a function you can plot and evaluate:

sol=First[y/.NDSolve[{y''[t]+y[t]^3==3Cos[t],y[0]==2,y'[0]==3},y,{t,0,100}]];
sol[t]

sol[1]

sol[101]

Plot[sol[t],{t,0,100}]

Plotting

Plot

plot1 = Plot[{Sin[x] + 2,Sin[x]},{x,-4,4},Frame->True,FrameLabel->{"x","y"},PlotStyle->Thick]

Filled Plots

plot1 = Plot[{Sin[x] + 2},{x,-4,4},Frame->True,FrameLabel->{"x","y"},PlotStyle->Thick,Filling->Axis]

plot1 = Plot[{Sin[x] + 2,Sin[x]},{x,-4,4},Frame->True,FrameLabel->{"x","y"},PlotStyle->Thick,Filling->{1}]

Combining Plots

plot1 = Plot[{Sin[x] + 2},{x,-4,4},Frame->True,FrameLabel->{"x","y"},PlotStyle->{Red,Dashed,Thick}]

plot2 = Plot[{Cos[x] + 2,Cos[x]},{x,-4,4},Frame->True,FrameLabel->{"x","y"},PlotStyle->Thick]

Show[plot1,plot2,PlotRange->{All, All}]

ParametericPlot

ParametricPlot[{(t^3),(t^2)/(1+(t^2))},{t,-3,3},Frame->True,FrameLabel->{"x","y"},AspectRatio->1,PlotStyle->Thick]

Plots of Implicit Functions

Mathematica plots implicit functions using the ContourPlot command:

ContourPlot[(x^4)-(y^6)==Sin[x*y],{x,-2,2},{y,-4,4},PlotPoints->100,ContourStyle->{Red,Thick}]

Plots of Explicit Families of Curves

You should pick values for the constants that provide a good picture of the family of curves. Expect to play around a bit with both the constants and the plot ranges until you get a picture that conveys the family of curves well.

For a family of curves with one constant:

y[t_]=Exp[-2t] + C[1]*((t^2)+1)

functionlist=y[t]/.C[1]->Range[-2,2,0.5]

lines=Plot[Evaluate[functionlist],{t,-2,2},PlotRange->{{-2,2},{-4,4}},AxesLabel->{"t","Family of Curves"},PlotStyle->Red]

For a doubly infinite family of curves:

y[t_]=C[1]Exp[-2t] + C[2]*Exp[-t]+t Sin[t]-1

functionlist=y[t]/.C[1]->Range[-2,2,1]/.C[2]->Range[-2,2,1]

lines=Plot[Evaluate[functionlist],{t,-2,6},PlotRange->{All,{-10,10}},AxesLabel->{"t","Family of Curves"},PlotStyle->Red]

Plots of Implicit Families of Curves

Probably the easiest way to plot a family of curves in Mathematica for an implicit function is to create a list of some members of the family, and then use ContourPlot to plot that list. The Evaluate command is sometimes needed to get the plot, so I included it here.

The Do command creates the list, so adjust the values of c to make sense for your problem (always choose some positive and some negative, since if you choose all positive you might miss interesting members of the family of curves). I have included output of sol so you can see the actual curves being sketched; if it is a long list, supress output using a semicolon. If the list is too long, it will take awhile for the plot to be created.

sol={};
Do[sol=Append[sol,x^2+y^3==c],{c,-2,2,0.5}]
sol

ContourPlot[Evaluate[sol],{x,-5,5},{y,-2,2},ImageSize->{400,Automatic},AspectRatio->0.5,FrameLabel->{"x","y"},PlotPoints->200]

Notice if you run your mouse over the curves, Mathematica will tell you the equation of each!

You can skip the table and let Mathematica choose the values of c for you, but frankly I prefer to choose my own c values (I feel the sketches are better). Notice that I turn off the contour shading when I do it this way. If you do it this way, Mathematica shows you the value of the contour when you mouse over, not the full equation.

ContourPlot[x^2+y^3,{x,-5,5},{y,-2,2},ImageSize->{400,Automatic},AspectRatio->0.5,FrameLabel->{"x","y"},PlotPoints->200,ContourShading->False]

Manipulate

Basic Example

Manipulate is a wonderful command that produces output as a function of a parameter. It can be used to create animations, among other things.

Plot[Cos[x],{x,0,2Pi},PlotStyle->Thick,PlotRange->{{0,2Pi},{-5,5}}]

We can investigate the behaviour of the function f(x)=cos(a x+b)+c using Manipulate.

Manipulate[Plot[Cos[a*x+b]+c,{x,0,2Pi},PlotStyle->Thick,PlotRange->{{0,2Pi},{-5,5}}],{a,-2,2},{b,0,2Pi},{c,-2,2}]

If you click on the plus sign beside each slider, you can see the value of the parameter that the plot is for. You can use the slider to change the value of the parameter, and even create animations! Try it! Obvioulsy, too many parameters still becomes a bit much if you are doing animations, but certainly this can help you explore the behaviour of fucntions relative to a parameter.

Example from Differential Equations

For example, here is a harmonic oscillator (starting at rest from the origin) as a function of mass (m) and spring constant (k) with a sinusoidal forcing function.

Clear[x,displacement]
sol = DSolve[{m x''[t]+k x[t]==Sin[t] ,x[0]==0,x'[0]== 0},x[t],t];
displacement[t_,k_,m_]=x[t]/.Last[sol]

We can easily  examine the baviour of the solution for  varying mass and spring constant.

Manipulate[Plot[displacement[t,k,m],{t,0,80},PlotRange->{{0,80},{-5,5}},PlotStyle->Thick,AxesLabel->{"t","Harmonic Oscillator"}],{m,1,5},{k,2,5}]

Notice that there are some value of k and m for which there is no plot. That is because the solution Mathematica found to the differential equation is not correct for those values of k and m. This is where k=m, and the solution Mathematica found is only valid where k≠m.

We can actually use Manipulate to invesitgate this further (notice I have now told Manipulate to vary the m and k by an integer amoun belwt). Notice that the form of the solution changes whenever k=m (you get a t factor in some of the terms, instead of pure trig solutions). Obvioulsy, we study differential equations to help us understand why this is happening.

Clear[x]
Manipulate[ DSolve[{m x''[t]+k x[t]==Sin[t] ,x[0]==0,x'[0]== 0},x[t],t],{m,1,5,1},{k,2,5,1}]

Example from Calculus III

The three types of torus--with the top cut off so you can see them.

Manipulate[ParametricPlot3D[{(c+1*Cos[t])Cos[s],(c+1*Cos[t])Sin[s],1*Sin[t]},{t,Pi,2Pi},{s,0,2Pi},PlotRange->{{-4,4},{-4,4},{-1,0}}],{c,0,2}]

Recursion

If you have used Mathematica before, you have probably encountered recursion. Recursion is easily set up in Mathematica, but you need to be careful or you can run into difficulties.

Let's say we want to numerically calculate a root of f(x)=2x cos(2x) - (((x-2))^2) that is near x=2 using Newtons method. We can do that recursively using the following:

f[x_]=2x Cos[2x]-(x-Pi)^2;
x[0]=2.0;
x[i_]:=x[i-1]-f[x[i-1]]/f'[x[i-1]]

data=Table[{i,x[i]},{i,0,10}];
TableForm[data]

This worked, but it is not a very good implementation of Newton's method. You should notice that it is somewhat slow, and that it only produced 5 decimals of output. We can get more output using SetPrecision , and speed it up by storing the intermediate results so Mathematica doesn't have to calculate them again and again.

f[x_]=2x Cos[2x]-(x-Pi)^2;
x[0]=2.0;
x[i_]:=x[i]=SetPrecision[x[i-1]-f[x[i-1]]/f'[x[i-1]],20]

data=Table[{i,x[i]},{i,0,10}];
TableForm[data]

There is one more important change to make when we are doing recursion--we want to clear the variables before each implementation. We can do that by adding a Clear command at the beginning of the cell. If we don't do this, it is possible that if you come back and re-evaluate the cell you will get errors due to some of the quantities being previously defined. Worse than getting errors, you can sometimes get no errors but get incorrect results--this will happen if you change the function f!

Clear[f,x,data]
f[x_]=2x Cos[2x]-(x-Pi)^2;
x[0]=2.0;
x[i_]:=x[i]=SetPrecision[x[i-1]-f[x[i-1]]/f'[x[i-1]],20]

data=Table[{i,x[i]},{i,0,10}];
TableForm[data]

Products and Sums

Mathematica can do some interesting things with products and sums. This is particularly useful when solving differential equations using infinite series, where you need to figure out the product form from the pattern, which can then be converted into a known function using Mathematica.

For example, let's say we have the function given by an infinte series with coefficients given by a set of recursion relations:

f(x)=(∑_(n=0)^∞)(Subscript[a, n])(x^n)
(Subscript[a, 0])=π
(Subscript[a, n])=((1×2×3×…×n)/(1×3×5×7×…×(2n-1)))(Subscript[a, 0])

Here's how we can figure out a simpler form for the function that does not involve recursion or summation. First, we use the Product command to determine (Subscript[a, n]):

a[n_]=Product[i/(2i-1),{i,1,n}]*Pi

Now, we use this value of (Subscript[a, n]) in the infinite sum for f(x):

f[x_]=Sum[a[n]x^n,{n,0,Infinity}]

This won't always work, but you will be surprised how often it does!

Built-in Function: Heaviside 
(used with Laplace Transforms)

Plot[HeavisideTheta[t],{t,-4,4},PlotStyle->{Red,Thick}]

Manipulate[Plot[HeavisideTheta[t-c],{t,-4,4},PlotStyle->{Red,Thick},PlotRange->{{-4,4},{0,1}}],{c,-4,4}]

Manipulate[Plot[HeavisideTheta[t-c]*Cos[t],{t,-4,4},PlotStyle->{Red,Thick},PlotRange->{{-4,4},{-1,1}}],{c,-4,4}]

Manipulate[Plot[(1-HeavisideTheta[t-c])*Cos[t],{t,-4,4},PlotStyle->{Red,Thick},PlotRange->{{-4,4},{-1,1}}],{c,-4,4}]

Manipulate[Plot[HeavisideTheta[t-c]*Cos[t-c],{t,-4,4},PlotStyle->{Red,Thick},PlotRange->{{-4,4},{-1,1}}],{c,-4,4}]

Built-in Function: Dirac delta 
(used with Laplace Transforms)

Plot[DiracDelta[t],{t,-4,4},PlotStyle->{Red,Thick}]

Clear[f]
Integrate[DiracDelta[t]f[t],{t,-20,20}]

Integrate[DiracDelta[t]f[t],{t,0,20}]

Integrate[DiracDelta[t-4]f[t],{t,0,20}]

Matrix Manipulations

Mathematica is quite comfortable working with matrices and vectors. The command MatrixForm is used to make the output look nice.

Here is a matrix:

A={{1,-1,2},{-1,-1/10,4},{6,3,0}}
MatrixForm[A]

And some vectors:

vec1={1,4,-3}
MatrixForm[vec1]

vec2={3,4,3}
MatrixForm[vec2]

Vector Addition

vec1+vec2

Dot Product

MatrixForm[A.vec1]

MatrixForm[Dot[A,vec1]]

Cross Product

MatrixForm[Cross[vec1,vec2]]

Matrix Power

MatrixPower[A,2]
MatrixForm[%]

Exponential of a Matrix

MatrixExp[N[A]]
MatrixForm[%]

Eigenvales and Eigenvectors

Eigenvalues[A]//N
Eigenvectors[A]//N
Eigensystem[A]//N

Complex Numbers

Mathematica is quite comfortable working with complex numbers.

Here are a couple of complex numbers:

Clear[a,b,x]

a=3-I
b=7+5I

The ComplexExpand function helps us break a complicated complex valued function into its real and imaginary parts.

ComplexExpand[Exp[a*x]]

ComplexExpand[Sin[b*y]Exp[a*x]]

Mathematica's programming syntax

We will not be doing a tremendous amount of programming in this class, but you will probably want to know how to some simple programming in Mathematica to help you with some of your assignments.

Looping

Do, While and For do not print any output; they simply do what they are asked to do. After the calculation you will need to output whatever you are interested in. Table is very similar to Do, but it forms a list of the output. In many situations Table is all you will need. Other commands that may prove useful are Sum and Product.

Do Loops

Clear[a,b,c]

a=0;
b=0;
c=0;
Do[
a=a+n;
b=a+b;
c=b^2,
{n,1,4}]

{a,b,c}

While Loops

Clear[a,b,c,n]

a=b=c=0; n=1;
While[n≤4,
a=a+n;
b=a+b;
c=b^2;
n=n+1]

{a,b,c}

For Loops

Clear[a,b,c,n]

For[a=b=c=0; n=1,
n≤4,
n=n+1,
a=a+n;
b=a+b;
c=b^2]

{a,b,c}

Table Loops

Clear[a,b,c]

a=b=c=0;
Table[
a=a+n;
b=a+b;
c=b^2,
{n,1,4}]

{a,b,c}

Branching

Branching in Mathematica can be accomplished via the commands If, Which, or Switch. The most common one you will use will probably be IF.

If

f[x_] = If[x≤0,x^2,x];
Plot[f[x],{x,-3,3},AxesLabel->{"x","y"},PlotStyle->{Blue}]

Flow Control

You can control the flow of programming in Mathematica using the commands Continue, Break, Return, and Goto. Here is a basic example which shows how these work.

f[n_] := 
Do[
Which[
n==1,Print["n=1 ","k= ",k]; Continue[],
n==2,Print["n=2 ","k= ",k];Break[],
n==3,Return["That's All Folks"],
True,Goto[jump];
Label[jump];Return[Print["n= ",n]]]
,{k,1,3}]

f[1]

f[2]

f[3]

f[4]

Selecting Elements from a List

To select elements from a list, use the command Part.

Clear[a,b,c]
data = Table[{a[n],b[n],a[n]*b[n]^2},{n,0,5}];
TableForm[data,TableHeadings->{None,{"1","2","3"}}]

data2 = Table[{Part[data,n+1,2],Part[data,n+1,3]},{n,0,5}];
TableForm[data2,TableHeadings->{None,{"2","3"}}]

A shorthand does exist for Part, and it is done by typing <esc> [[ <esc> and then the part of the item you want followed by <esc> ]] <esc>. Here is an example, try it yourself!

The full list of data:

data

The third element in data:

data[[3]]

The second element in the third element in data:

data[[3,2]]

Numerical Output

If you simply leave Mathematica to its defaults, it will output numbers to 6 digit precision, or significant digits precision. For example:

N[Pi]
N[Pi,10]
N[3/4]
N[3/4,8]

If we want to include more digits in the output, we should use the command SetPrecision when we output the number. We would still like Mathematica to work with high digits of precision internally, but we want to be able to see as many digits as we like when needed.

SetPrecision[Pi,9]
SetPrecision[3/4,1]
SetPrecision[3/4,5]

The machine precision of Mathematica can be determined:

$MachinePrecision

If we use the symbolic representation of numbers in Mathematica to do arithmetic (variable precision arithmetic) we get no round off error, but our calculations will be slow.

Pi-Pi
Timing[Sum[(((1-24)/n))^n ,{n,1,600}]];
SetPrecision[%,16]

That took over 2 seconds on my computer (the first number), and the result was -1.3189...(we actually calculated it as a rational number, not a decimal, but it is so long I supressed the output!).

If we use the hardware to do our arithmetic (fixed precision arithmetic), we get results much faster, but we introduce roundoff error:

N[Pi,25]-N[Pi,25]
Timing[Sum[(((1-24.0)/n))^n ,{n,1,600}]]
SetPrecision[%,16]

This may not seem like a big time difference for the accuracy you get out at the end, and it really isn't for this simple case. However, when you inadvertantly use arbitrary precision numbers inside a looping structure, you can end up with code which takes days to run instead of minutes. The lesson: know your expected error and work within it!

Putting It All Together: Euler's Method

For solutions to assignment problems, I am not expecting you to do high level programming in Mathematica. For example, here are some ways of implementing Euler's method to solve an initial value problem. I would expect you to be able to create the recursion relation method, or the While loop method, but not the module method. However, you should feel free to explore some of the more advanced syntax of Mathematica if it pleases you!

Here we want to solve the initial value problem y'(t) = f(y,t), y(α)=β for y(t) in a given range β≤t≤γ. Euler's method generates approximations to the solution y(t) using the following relations, once you have picked a value for the step size h:

(Subscript[t, n])=(Subscript[t, n-1])+h
(Subscript[y, n])=(Subscript[y, n-1])+h *f((Subscript[y, n-1]),(Subscript[t, n-1])), n=1,2,3,...
(Subscript[t, 0])=α,(Subscript[y, 0])=y(α)=β.

Let f(y,t)=y-(t^2)+1, y(0)=1/2 and 0≤t≤2. I've included some extra tables just to show you what Mathematica can do.

Euler's Method: While Loop

Clear[f,y,t,a,b,α,data,num]
f[t_,y_]=y-(t^2)+1;
a=0.0;
b=2.0;
α=0.5;
num=10;
h=((b-a)/num);
yi=α;
ti=a;
data={{ti,yi}};
n=1;
While[n≤num,{ti,yi}=N[{t+h,y+h f[t,y]}/.{t->ti,y->yi}];data=Append[data,{ti,yi}];n=n+1]
SetPrecision[data,9]
ListPlot[data,PlotRange->{{0,2},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Blue,Dashed}]

Euler's Method: Recursion Relations

Clear[f,y,t,a,b,α,data]
f[t_,y_]=y-(t^2)+1;
a=0.0;
b=2.0;
α=0.5`;
num=10;
h=((b-a)/num);
t[0]=a;
y[0]=α;
y[n_]:=y[n]=y[n-1]+h f[t[n-1],y[n-1]]
t[n_]:=t[n]=a+n h
data=Table[{n,t[n],N[y[n],20]},{n,0,num}];
SetPrecision[data,9]
data2=Table[{data[[n+1,2]],data[[n+1,3]]},{n,0,num}];
SetPrecision[data2,9]
phase=Table[{data[[n+1,2]],data[[n+1,3]]},{n,0,num}];
ListPlot[phase,PlotRange->{{0,2},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Red,DotDashed}]

Euler's Method: Modules

Modules allow you to create executable blocks of code, with variables that are internal to the module.

Clear[f,y,t,a,b,α,data]

euler[f_,t_,y_,t0_,y0_,h_,n_]:=Module[{ti=t0,yi=y0},Prepend[Table[{ti,yi}=N[{t+h,y+h*f} /.{t->ti,y->yi}],{k,1,n}],{t0,y0}]]

Here is the output from executing the module:

euler[ y-t^2+1,t,y,0.0,0.5,0.2,10]

Here the output is put into a table, and plotted.

approx=euler[y-(t^2)+1,t,y,0.,0.5,0.2,10];
SetPrecision[approx,8]
ListPlot[approx,PlotRange->{{0,2},{0,5}},Joined->True,Frame->True,FrameLabel->{"t","y"},PlotStyle->{Dotted}]