8.Maths - Mathematics for Physicists

Introduction

In these examples, we have put together a very incomplete crash course in Mathematica and try to give you a flavour for some of the things you can do with Mathematica. We also try to emphasise ways in which it can go wrong or be misused. Always think about the results you get from a computer and try to find ways to verify them, for example by thinking about limiting cases. The Mathematica help system is indeed very helpful and gives you a lot of information and examples (some of which are used here) about the functions we  discuss.

0.1 Functions and equations

Mathematica provides many in-built functions (always with a capitalized first letter) that can be either symbolically or numerically manipulated

Sin[x]

Sin[5.3]

Many identities are built in and implemented through the Simplify[] or FullSimplify[] or TrigExpand[] functions

TrigExpand[Sin[2x]]

Simplify[Sin[2x]-2Sin[x]Cos[x]]

Many polynomials and special functions are defined

LegendreP[5,x]

WhittakerM[4,3,5.4]

Selecting a function name and pressing F1 brings up the inbuilt help on that function

You can also define your own functions  ( the  := is a delayed assignment)

f[x_]:=(3/(4x+Sin[x]))

f[4.5]

f[y]

Mathematica can solve equations or sets of equations in various different contexts

Solve[2(x^3)+(x^2)+a x+1==0,x]

Solve[{A x + y == 7 , B x - y == 1}, {x, y}]

When an analytic solution is not available, Mathematica can use numerical methods

NSolve[2(x^3)+(x^2)+7 x+1==0,x]

Solve the equation sin(z+sin(z+sin(z)))==cos(z+cos(z+cos(z))):

roots=NSolve[Sin[z+Sin[z+Sin[z]]]==Cos[z+Cos[z+Cos[z]]]&&-3<Re[z]<3&&-3<Im[z]<3,z];
ListPlot[{Re[z],Im[z]}/.roots,PlotLabel->TraditionalForm[Sin[z+Sin[z+Sin[z]]]==Cos[z+Cos[z+Cos[z]]]]]

0.2 Derivatives and Integrals

Derivatives are taken using D[f,x] or D[f,{x,n}]

D[Sin[x],x]

Higher order derivatives

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

Partial derivatives too (((∂^2)f)/(∂x∂y))

D[Sin[x y], x,y]

Indefinite and definite integrals (don’t forget integration constants)

Integrate[Sin[x],x]

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

There are many built in integration rules that Mathematica can use (almost all of what you find in tables of integrals is in Mathematica)

Integrals can also be evaluated numerically

NIntegrate[(x^(x Log[x-2])),{x,-2,5}]

Mathematica uses some sophisticated techniques to integrate. It  automatically chooses what method it uses but the user can specify the method too - this is useful but can also lead to no being able to integrate precisely. In this latter case, Mathematica indicates it is unable integrate

Plot[(1/x)Cos[Log[x]/x],{x,0,1},PlotRange->15]
NIntegrate[(1/x)Cos[Log[x]/x],{x,0,1}]
NIntegrate[(1/x)Cos[Log[x]/x],{x,0,1},Method->"GlobalAdaptive"]
NIntegrate[(1/x)Cos[Log[x]/x],{x,0,1},Method->"AdaptiveMonteCarlo"]
NIntegrate[(1/x)Cos[Log[x]/x],{x,0,1},Method->"NewtonCotesRule"]

0.3 Replacement rules and list manipulations

It is often useful to make replacements in expression in Mathematica

The % symbol indicates the previous output and %n indicates output number n

Expressions can be easily expanded (Expand[]) and less easily simplified (Simplify[] or FullSimplify[])

3(x^3)+2x+9 /. x->(x-3)
Expand[%]
FullSimplify[%]

Replacements can be repeated using the //. notation (this can lead to recursive strife)

3(x^3)+2x+9 //. x->(y-3)

3(x^3)+2x+9 //. x->(x-3)

A fundamental object in Mathematica is a list

list={3,4,a,7-2x,2(2+a x)}

Lists can be accessed element wise or as sublists

list[[2]]

list[[2;;4]]

Lists of lists

listsoflists={{31,14,(a^2)},{17-2x,5y},2(2+a x)}

Length[listsoflists]

listsoflists[[1,3]]

Functions can be applied to lists and act elementwise

Sin[list]
Sin[listsoflists]

0.4 Plotting

Plotting functions is easy in Mathematica and by specifying options carefully can make high quality graphs.

Plot[Sin[x],{x,0,2π}]

You can plot data

dat=Table[Around[RandomReal[5],RandomReal[0.5]],20]
ListPlot[dat]

You can also plot in 3D

Plot3D[(x^2 + y^2) Exp[-(x^2 + y^2)], {x, -2, 2}, {y, -2, 2}]

There are also many other ways to plot calculations or data

ContourPlot[(Sqrt[(((y z + 10))^2)+(y^2)+(z^2) ]),{y,-5,5},{z,-5,5}, ContourLabels->True]

1.1 Taylor expansions

The Series[] function is useful in this regard.

Series[f[x],{x,(Subscript[x, 0]),n}] provides the expansion of f around the point (Subscript[x, 0]) to order n

Series[] can construct standard Taylor series, as well as certain expansions involving negative powers, fractional powers, and logarithms.

ser=Series[(1/(1-x)),{x,0,10}]

ser2=Series[(1/(1-x-(x^2))),{x,0,10}]

x0=2;
Table[(1/(n!))D[(1/(1-x)),{x,n}]/.x->x0,{n,10}]

CoefficientList[ser2, x]

SeriesCoefficient[ser,10]

serxy=Series[Sin[x+y],{x,0,3},{y,0,3}]

Asymptotic behaviour can also be extracted using Limit

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

Limit[((1+Sinh[x])/(Exp[x])),x->∞]

Limit[Abs[x]/x, x -> 0, Direction -> "FromAbove"]
Limit[Abs[x]/x, x -> 0, Direction -> "FromBelow"]

Series[] outputs can be further manipulated which is useful for approximately solving

Expand[((Normal[ser])^2)]

Mathematica can also sum a lot of known infinite series

Sum[(n+1) (x^n),{n,0,Infinity}]

Sum[(1/(n^2)),{n,1,∞}]

1.2 First order ordinary differential equations

Mathematica can be used to solve differential equations either analytically or numerically

Here are a couple of examples you can solve by separation of variables. With only the DE specified, the general solution is given

eq1sol=DSolve[y'[x]==((3(x^2))/(y[x])),y,x]

eq2sol=DSolve[y'[x]==((Cos[x])/(Sin[2y[x]])),y,x]

If you want specific solutions you need to specify the boundary conditions too

eq3sol=DSolve[{y'[x]==((3(x^2))/(y[x])),y[3]==0},y,x]

You must specify the right number of boundary conditions!

DSolve[{y'[x]==((3(x^2))/(y[x])),y[3]==0, y'[2]==1},y,x]

We can easily plot such solutions (both are shown)

Plot[Re[y[x]/.eq3sol],{x,0,10}]
Plot[Im[y[x]/.eq3sol],{x,0,10}]

Mathematica has many ways of solving DEs, but there are plenty of cases where is fails to solve them (in which case it returns your input).

DSolve[y'[x]==Exp[((Cos[x])/(Sin[2y[x]]))],y,x]

Sometimes it is helpful to have a numerical solution

NDSolve[{eq,bcs},y,{x,xmin,xmax}] provides an approximate solution over the range {xmin,xmax}.

eq1nsol=NDSolve[{y'[x]==3((y[x])^2),y[1]==1},y,{x,0,1}]
Plot[{y[x]/.eq1sol,y[x]/.eq1nsol},{x,0,1}]

This is a case where an analytic solution is easily found

eq1sol=DSolve[{y'[x]==3((y[x])^2),y[1]==1},y,x]

The numerical solution agrees by eye in thee region where it was approximated, but is incorrect outside that region.

Plot[{y[x]/.eq1sol,y[x]/.eq1nsol},{x,-1,1}]

eq1nsol=NDSolve[{y'[x]==Exp[((Cos[x])/(Sin[2y[x]]))],y[1]==1},y,{x,0,1}];
Plot[y[x]/.eq1nsol,{x,0,1}]

eq1nsol=NDSolve[{y'[x]==Exp[((Cos[x])/(Sin[2y[x]]))],y[1.000000001]==1},y,{x,0,1}];
Plot[y[x]/.eq1nsol,{x,0,1}]

eq1nsol=NDSolve[{y'[x]==((3(x^2))/(y[x])),y[3]==0},y,{x,-5,5}]

NDSolve[{y'[x]==Exp[((Cos[x])/(Sin[2y[x]]))],y[1]==1},y,{x,-5,5}]

1.3 Second order ordinary differential equations

The same DSolve and NDSolve operations also work for second order equations

The default NDSolve finds a solution to the following DE

eqsol = NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0},x,t]
Plot[x[t]/.eqsol,{t,0,10}]

Internally NDSolve is automatically selecting one of many methods that it knows. You can specify the method to find solutions to  the same equation.

Experimenting with solution methods and seeing consistency can give you some confidence in the solution, although without uniqueness guarantees that confidence should be tempered.

Here is  case where many solutions are exist and we try various initial guesses for the slope at x=0.

eqsol =Table[NDSolve[{x''[t] + Sin[x[t]] == 0, x[0] == x[10] == 0}, x, t, Method -> "BoundaryValues" -> {"Shooting", "StartingInitialConditions" -> {x[0] == 0, x'[0] == dd}}],{dd,{1.5,1.75,2,-1.75}}]
Plot[x[t] /. eqsol, {t, 0, 10}]

1.4 General linear ordinary differential equations

The same tools work with more general ODEs

2.1 Two variables

Matrices in Mathematica are no more than lists of lists of the appropriate dimensions

vec={a,b,c,d}
vec//MatrixForm

mat={{3,4,7 I,6},{-2,4,(11/2),4},{22,-22/10,22,-2},{0,1,0,0}}
mat//MatrixForm

Note “I” is the imaginary number

Transpose[mat]//MatrixForm
ConjugateTranspose[mat]//MatrixForm

Some special matrices are built-in

IdentityMatrix[3]//MatrixForm

Table[PauliMatrix[i]//MatrixForm,{i,3}]

DiagonalMatrix[{2,5,I,h}]//MatrixForm

Matrix operations are pretty much fully supported. Matrix-vector or matrix matrix multiplication is through the . operator

mat.vec//MatrixForm

mat.mat // MatrixForm

Note that functions applied to matrices are component wise

(mat^2)//MatrixForm

There are built in routines for determinants

Det[mat]

For large numerical matrices, the determinant is clearly not computed as you would have done it in high-school.

n=1000;
Det[Table[RandomReal[],{i,n},{j,n}]]

you can use this to find eigenvalues

characteristiceqn=Det[mat -λ IdentityMatrix[4]]

Given the complex natrure of the solutions, finding then by inspection is not wasy

Plot3D[{Re[characteristiceqn/. λ->x+I y],Im[characteristiceqn/. λ->x+I y]},{x,-5,5},{y,-10,10}]

Solve[N[characteristiceqn]==0,λ]

There are also built in functions to do this (note the different ordering of the eigenvalues (

Eigenvalues[mat]
Eigenvectors[mat]
Eigensystem[mat]

These work with symbolic matrices too

anothervec={1,2,3,4};

(anothermat={{3a+b,-y,s},{3,a/b,y^2+a},{2,a,a}})//MatrixForm

Det[anothermat]

Eigenvalues[anothermat]

Eigenvectors can be used to diagonalise a matrix  (the N[] function gives the approximate numerical representation

evecs=Eigenvectors[N[mat]]

tvec=Transpose[evecs];
invtvec=Inverse[tvec];

(diag=invtvec.mat.tvec)//MatrixForm

This doesnt look diagonal, but it really is with machine rounding errors on the off diagonal. We can throw away those with Chop

Chop[diag]//MatrixForm

We can also redo the calculation in high precision (the N[] function specifies the precision with which to work)

evecs = Eigenvectors[N[mat,30]]
tvec=Transpose[evecs];
invtvec=Inverse[tvec];
(diag=invtvec.mat.tvec)//MatrixForm

2.2 Multiple variables

Functions of matrices must be dealt with carefully. Eg, the following gives the component-wise function

Exp[mat]//MatrixForm

MatrixExp[mat]//MatrixForm

MatrixFunction[Cos,mat]//MatrixForm

Lets solve the coupled set of ODEs ((dx)/(d t))=M x for a set of four functions

fvec={a[t],b[t],c[t],d[t]}

sol=DSolve[{D[fvec,t]==mat.fvec,{1,1,1,1}==(fvec/.t->0)},fvec,t]

The solution is exactly exp(M t) x(0) as

Chop[MatrixExp[t mat].{1,1,1,1}]

2.3 Higher order coupled linear ODEs

2.4 Symmetries in matrices

n=10;

S1=Table[0,{i,n},{j,n}];
Do[S1[[i,i+1]]=1,{i,n-1}];
S1[[n,1]]=1;
S1//MatrixForm

Chop[Eigenvectors[N[S1]]]

Tnprime=Table[0,{i,n},{j,n}];
Do[Tnprime[[i,i]]=2;
Tnprime[[i,i+1]]=-1;
Tnprime[[i,i-1]]=-1;
,{i,2,n-1}];
Tnprime[[1,1]]=Tnprime[[n,n]]=2;
Tnprime[[1,2]]=Tnprime[[1,n]]=Tnprime[[n,1]]=Tnprime[[n,n-1]]=-1;
Tnprime//MatrixForm;

evecs=Eigenvectors[Tnprime];

ListPlot[evecs[[90;;100;;1]],Joined->True,PlotRange->All]

Operations on matrices are elementwise by default

m = {{a, b}, {c, d}};
m^2

m.m

MatrixPower[m,2]

3.1 From particles to fields

3.2 Solving PDEs

General solutions to PDEs are given in terms of arbitrary functions (Subscript[c, i])

eqn = x^2*D[u[x, y], x] + y^2*D[u[x, y], y] - (x + y)*u[x, y];
sol = DSolve[eqn == 0, u, {x, y}]

Model the vibrations of a string with fixed length, say π, using the wave equation:

eqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}];

Specify that the ends of the string remain fixed during the vibrations:

bc={u[0,t]==0,u[π,t]==0};

Obtain the fundamental and higher harmonic modes of oscillation:

ic={u[x, 0] == Sin[m x], Derivative[0, 1][u][x, 0] == 0};

dsol=DSolve[{eqn,bc,ic},u,{x,t},Assumptions->Element[m,Integers]]

Visualize the vibrations of the string for these modes:

sol1[x_,t_]=u[x,t]/.dsol[[1]]

Table[Show[Plot[Table[sol1[x,t],{t,0,4}]//Evaluate,{x,0,Pi},Ticks->False]],{m,4}]

In general, the solution is composed of an infinite number of harmonics:

ic={u[x, 0] == x(π-x)^2, Derivative[0, 1][u][x, 0] == 0};

dsol=DSolve[{eqn,bc,ic},u,{x,t}]/.{K[1]-> m}

Extract four terms from the (RefLink[Inactive, paclet:ref/Inactive]) sum:

Clear[sol]

sol[x_,t_]=u[x,t]/.dsol[[1]]/.{∞->4}//Activate

Visualize the vibration of the string:

Animate[Plot[sol[x,t],{x,0,π},PerformanceGoal->"Quality",PlotRange->{-5,5},ImageSize->Medium],{t,0,12},SaveDefinitions->True]

Model the oscillations of a circular membrane of radius 1 using the wave equation in 2D:

eqn=r D[u[r,t],{t,2}]==D[r D[u[r,t],r],r];

Specify that the boundary of the membrane remain fixed:

bc=u[1,t]==0;

Initial condition for the problem:

ic={u[r,0]==0,Derivative[0,1][u][r,0]==1};

Obtain a solution in terms of Bessel functions:

(dsol=DSolve[{eqn,bc,ic},u[r,t],{r,t}])//TraditionalForm

Extract a finite number of terms from the (RefLink[Inactive, paclet:ref/Inactive]) sum:

h[r_, t_]=u[r,t]/.dsol[[1]]/.{∞->3}//Activate//N;

Visualize the oscillations of the membrane:

Animate[Plot3D[h[r,t]/.{r->Sqrt[x^2+y^2]},{x,y}∈Disk[],PlotRange->{-1,1},Ticks->None,Mesh->True,MeshStyle->{Red, Blue},PlotStyle->Yellow],{t,0,4},SaveDefinitions->True]

When analytic methods are not know (eg for strange geometries), numerical solutions can be helpful

Simple model for soil temperature at depth x with periodic heating at the surface:

NDSolve[{D[u[t,x],t]==D[u[t,x],x,x],u[0,x]==0,u[t,0]==Sin[t],u[t,5]==0},u,{t,0,10},{x,0,5}]

Plot3D[Evaluate[u[t,x]/.%],{t,0,10},{x,0,5}, PlotRange->All]

Wave equation for geometry with a slit (might take a while to run!)

Ω=RegionDifference[RegionDifference[Rectangle[{0,0},{2,1}],Rectangle[{9/10,0},{11/10,4/10}]],Rectangle[{9/10,6/10},{11/10,1}]];
sol=u/.First[NDSolve[{D[u[t,x,y],{t,2}]-Laplacian[u[t,x,y],{x,y}]==0,DirichletCondition[u[t,x,y]==0,True],u[0,x,y]==2*Exp[-125((x-0.25)^2+(y-0.5)^2)],Derivative[1,0,0][u][0,x,y]==0},u,{t,0,2},{x,y}∈Ω]]

ListAnimate[Table[Rasterize[Plot3D[sol[t,x,y],{x,y}∈Ω,PlotRange->{-0.75,2},AspectRatio->Automatic]],{t,0,2,1/25}],SaveDefinitions->True]

3.3 Fourier analysis

Fourier series can be computed, but are also built in

Here is a function

f1[x_]=UnitStep[Mod[x+π/2,2π,-π]]UnitStep[Mod[π/2-x,2π,-π]]

Plot[f1[x],{x,-3π,3π},Exclusions->None]

Now lets approximate it with a Fourier series

FourierSeries[f1[x],x,6]

We can compute the coefficients directly

coeff[m_]=FullSimplify[Integrate[f1[x]((Exp[-I m x])/(2π)),{x,0,2π}]]

For m=0, this will complain about 0 x ComplexInfinity, so here is a workaround

zeroval=Integrate[((f1[x])/(2π)),{x,0,2π}]

The “_” in this denotes the default (so if m=0 this returns zerovval  and otherwise it returns coeff[m]

carefulcoeff[m_]:=Switch[m,0,zeroval,_,coeff[m]]

These reproduce the Fourier series function above as we would hope

Table[{m,carefulcoeff[m]},{m,-5,5}]

Now we can define successively better approximations to the original function

fsapprox[Nmax_,x_]:=Sum[Exp[I m x]carefulcoeff[m],{m,-Nmax,Nmax}]

fsapprox[3, x]

Plot[{f1[x],fsapprox[3,x],fsapprox[10,x],fsapprox[100,x]}, {x, -3 π, 3 π}, Exclusions -> None]

We can make a nice animation of this

Manipulate[Plot[{f1[x],fsapprox[m,x]}, {x, -3 π, 3 π}, Exclusions -> None],{m,0,100,1}]

Fourier transforms are also implemented for many types of functions (as are other transforms such as Laplace and Hankel transforms)

fsaw[t_]=Piecewise[{{-2(a+t),-a≤ t<0},{2(a-t),0<t<a}}];

Plot[fsaw[t]/.a->1,{t,-2,2}]

ftsaw[α_]=FullSimplify[FourierTransform[fsaw[t],t,α,FourierParameters->{-1,-1}],Assumptions->{a>0}]

iftftsaw[t_]=FullSimplify[InverseFourierTransform[ftsaw[α],  α,t, FourierParameters -> {-1, -1}]]

Plot[fsaw[t] /. a -> 1, {t, -2, 2}]

GraphicsGrid[{{Plot[fsaw[t]/.a->1,{t,-2,2}],Plot[Im[ftsaw[α]/.a->1],{α,-15,15}],Plot[iftftsaw[t]/.a->1,{t,-2,2}]}}]

3.4 Scalar fields in higher dimensions

3.5 Vector fields