Mathematica: Getting Started
Branislav K. Nikolić
Department of Physics & Astronomy, University of Delaware, Newark, DE 19716
EMAIL: bnikolic@udel.edu

1. Basic Conventions

ClearAll["Global`*"] (*erase all variables*)

Cos[2] (*One of the first stumbling points in Mathematica is the use of parentheses and brackets.In Mathematica, brackets[] are always used to indicate the argument of a function*)

Cos[2]//N (*produce numerical value*)

cos[2]//N  (*built in Mathematica functions are always capitalized. Forgetting to capitalize Cos produces no result*)

N[Cos[2],1000] (*numerical computation can be done with arbitrary precision*)

Parentheses ( ) cannot be used to delineate arguments of  a function, instead use them to group arithmetic expressions:

N[8*Pi/(2 + 3),10]

(4-a) * 6

1.2 Defining functions
When defining a function you are defining an operation that can be applied to different objects in Mathematica. In defining a function, you have to remember two basic syntax rules: you have “_” after each variable in brackets, and you have “:” before the equals sign

F[x_]:=x^2

F[5]

g[r_,s_]:=r^2*(s-2)

g[F[3],F[5]]

1.3 Lists
To be able to easily manipulate a large number of objects at once, using a list is necessary. {} curly braces are used in Mathematica to denote a list:

x={3,5,7}

y={a,b,c}

After you have defined a list, it is easy to carry out operations on all of its elements at once.

x+2

x-y

1.4 Equality operators = vs. ==
Most of the equality statements you use in Mathematica will be simple definitions like:

x=3

y=2+x

Remember when defining variables that Mathematica is case-sensitive; there's a difference between vel and Vel. However, sometimes rather than just defining the equality, you will want to test two expressions to find out if they are equal. You can do this in Mathematica by using two equal signs:

1==1

1==2

A more frequent use of the double equals sign is in the Solve expressions, where the computer uses the equality test to solve equations:

Solve[2*p == 6, p]

Forgetting the second equals sign in the Solve expression is a common mistake which is easy to miss when you are searching for your error.

1.5 % operator
If you want to use the last result you obtained in the next expression, you can use % to represent it:

2+3

2*%

This can be a convenient shorthand. However, it is important to remember if you are debugging and rerunning parts of your program, that the % refers to the last result Mathematica obtained, and not necessarily to the result appearing above the next expression on the screen.

1.6 Substitution rule
The substitution rule is a very powerful tool in Mathematica. It removes the need to constantly redefine variables to substitute into expressions. A substitution rule is written like this:

	Expression /. Variable -> value

x=.  (*clear any value of variable*)

x

x^2/.x->2

x

By doing this you performed the operation without permanently assigning 2 to x. Another valuable use for substitution rules is in working with solutions to equations:

sol = Solve[x^2-3*y == 0, x]

The solution to this equation is contained in a nested list, which is a list of lists. Lists in the Wolfram Language are represented by {}. Items in lists are called elements and can be referred to by their position.To use the solution, you must first get it out of the nested list. Use [[ ]] (the short form for Part), with the position of the solution in the nested list. The solution here is the first (and only) element in the nested list:

sol[[1]]

x^2-3*y==0/.sol[[1]] (*plug only first solution into original equation*)

x^2-3*y==0/.sol (*plug both solutions into original equation*)

x^2/.sol

x^2

Clear["Global`*"] (*erase all previous variables*)

eqns={x+3==4,y-5==2,3 z==4} (*solve system of equations*)

u=Solve[eqns,{x,y,z}] (*solution comes out as a list*)

u[[1]] (*remove outer curly brace*)

{x->1,y->7,z->(4/3)}

eqns/.u[[1]]  (*check if equations are properly solved*)

x+y+z/. u[[1]] (*plugs in all three solutions into an expression*)

u[[1,1]] (*select only the first solution as substition rule for future usage*)

u[[1,2]]

x+y+z/. u[[1,1]] (*plug in only the first solutions*)

2. Numbers

In 1903 at a meeting of the American Mathematical Society, F. N. Cole read a paper entitled “On the Factorization of Large Numbers.” When called upon to speak, Cole walked to the board and, saying nothing, raised two to its sixty-seventh power and subtracted one from the answer. Then he multiplied, longhand, 193,707,721 by 761,838,257,287 and the answers agreed. Without having said a word, Cole sat down to a standing ovation. Afterwards he announced that it had taken him twenty years of Sunday afternoons to factorize the Mersenne number (2^67)-1.

2^67-1

FactorInteger[2^67-1]//Timing

FactorInteger[2^227-1]//Timing

?FactorInteger

??FactorInteger

?Prime*

1260/1530//N

N[1260/1530]

N[Pi,1000]

N[E,100]

1/Infinity

a=Binomial[6,2]

b=(1/6)^2*(5/6)^4

c=a*b

a

a=. (*clear any assigned value of a specific variable*)

a^2

Clear[a,b,c] (*clear any specific value of several variables*)

b

3. Elementary Functions

exponential                              Exp[x]
natural logarithm                  Log[x]
logarithm to base b            Log[b,x]
trigonometric                          Sin, Cos, Tan, Sec, Csc, Cot, ArcSin, ArcCos, ArcTan, ArcSec, ArcCsc, ArcCot
inverse trigonometric          ArcSin, ArcCos, ArcTan, ArcSec, ArcCsc, ArcCot
hyperbolic trigonometric   Sinh, Cosh, Tanh, Sech, Csch, Coth
inverse hyperbolic trig         ArcSinh, ArcCosh, ArcTanh, ArcSech, ArcCsch, ArcCoth

Log[10,1000]

Log[E]

Log[E] (*type Euler number as ESCAPE ee ESCAPE*)

Sin[Pi/6]

Sin[π/6] (*Type Green letter for Pi using ESCAPE pi ESCAPE*)

Sin[30 Degree]

ArcCos[1/Sqrt[2]]

expr=Sin[x]*Cos[x]

expr/.x->Pi/5//N

a=Table[Random[Integer,{1,6}],{10}]

Max[a]

4. Simplifying Expressions

4.1 Simplifying algebraic expressions

a=(1+x)*(2+x)+(1+x)^2

a//Simplify

a//Factor

(1+x) (3+2 x)

a//Expand

n!/(n-1)!//Simplify

n!/(n-1)!//FullSimplify

(1+x+x^2-x^3)/(x+2)^2//Apart

?Apart

4.2 Simplifying nonalgebraic expressions

Simplify[Gamma[x] Gamma[1-x]] (*No progress with Simplify*)

FullSimplify[Gamma[x] Gamma[1-x]](*FullSimplify tries a much wider range of transformations,involving not only algebraic functions,but also many other kinds of functions*)

4.3 Simplifying with Assumptions

Simplify[2*a+2*Sqrt[a-Sqrt[-b]]*Sqrt[a+Sqrt[-b]]]

FullSimplify[2*a+2*Sqrt[a-Sqrt[-b]]*Sqrt[a+Sqrt[-b]]] (*Im this case, brute force usage of FullSimplify does not bring any progress*)

Simplify[2*a+2*Sqrt[a-Sqrt[-b]]*Sqrt[a+Sqrt[-b]],Assumptions->{a>0&&b>0}] (*Key trick is to add Assumptions*)

Simplify[ArcSin[Sin[x]]]

FullSimplify[ArcSin[Sin[x]]](*Im this case, brute force usage of FullSimplify does not bring any progress*)

Simplify[ArcSin[Sin[x]],Assumptions->{-Pi/2<x<Pi/2}] (*Key trick, again, is to add Assumptions*)

4.4 Simplifying trigonometric expressions

TrigExpand — Sin[2x]⟶2 Sin[x]Cos[x] etc.
TrigFactor — Sin[..]+Sin[..]⟶product of terms, etc. — factors a trigonometric expression into a product of terms
TrigReduce — 2 Sin[x]Cos[x]⟶Sin[2x] etc. — rewrites products and powers of trigonometric functions in terms of trigonometric functions with combined arguments
TrigToExp, ExpToTrig — convert between trigonometric functions and exponentials
ComplexExpand — expand into real and imaginary parts

?Trig*

Sin[ArcSin[t]+ArcSin[v]]//TrigExpand

RegionPlot[ArcSin[Sqrt[1-t^2] v+t Sqrt[1-v^2]]==ArcSin[t]+ArcSin[v],{t,-1,1},{v,-1,1}] (*because of the multivaluedness of ArcSin this identity does not hold everywhere*)

TrigFactor[Sin[x]^2+Tan[x]^2]

Simplify[(1+Sin[x]-Cos[x])/(1+Sin[x]+Cos[x])]

Simplify[Sin[x]^3*Cos[x]^3] (*will not produce anything*)

FullSimplify[Sin[x]^3*Cos[x]^3] (*will not produce anything*)

TrigReduce[((Cos[x])^3) ((Sin[x])^3)] (*using dedicated TrigReduce will simplify expression by rewriting products and powers of trigonometric functions in terms of trigonometric functions with combined arguments*)

5. Solving algebraic equations

5.1 Solving algebraic equations symbolically

Clear[a,b,c,x,y]

eqn=a* (x^2)+b *x+c==0

sol=Solve[eqn,x]

x/.sol

eqn/.sol//Simplify

sol[[1]]

eqn/.sol[[1]]//Simplify

Reduce[a (x^2)+b x+c==0,x]

?Solve

?Reduce

Solve[{2x+5y==4,x-3y==3},{x,y}]

Solve[Cos[x]==0,x]

Reduce[Cos[x]==0,x]

Solve[x^2+y^2≤r^2,{x,y},MaxExtraConditions->All]

Reduce[x^2+y^2≤r^2,{x,y}]

5.2 Solving algebraic equations numerically

eqn2=a*x^4-b*x^3+c*x^2-1==0

Solve[eqn2,x]

eqn3=4*x^4-3*x^3+2*x^2-1==0 (*switch to numerical, instead of symbolic, coefficients, as general formulas are too cumbersome*)

Solve[eqn3,x]

sol1=NSolve[eqn3,x]

sol1[[3]]

Solve[Exp[a*x]==b,x]

f=Exp[-x]-0.5*x

Plot[f,{x,0,3}]

N[FindRoot[f==0,{x,0.5}],100]

f/.%

6. Plots of Real and Complex Functions

6.1 2D Plots of Real Functions

Plot[Exp[-x]*Sin[2x],{x,0,2π}]

Plot[{Sin[x],Cos[x]},{x,0,2π},PlotLegends->"Expressions"]

6.2 3D Plots of Real Functions

f=Sin[x^2]*Cos[Sqrt[y]]

Plot3D[f,{x,0,3},{y,0,4}]

ContourPlot[f,{x,0,3},{y,0,4},PlotPoints->40]

ContourPlot3D[Sin[x y+10^16]==z,{x,-2,2},{y,-2,2},{z,-1,1}] (*Machine-precision arithmetic*)

ContourPlot3D[Sin[x y+10^16]==z,{x,-2,2},{y,-2,2},{z,-1,1},WorkingPrecision->20] (*Arbitrary-precision arithmetic*)

ContourPlot3D[x^4+y^4+z^4-(x^2+y^2+z^2)^2+3 (x^2+y^2+z^2)==3,{x,-2,2},{y,-2,2},{z,-2,2}] (*rotate yourself to emphasize features*)

SphericalPlot3D[Sin[θ],{θ,0,Pi},{ϕ,0,2 Pi}]

RevolutionPlot3D[x^4-x^2,{x,-1,1}]

6.3 Plots of Complex Functions

Plot3D[Im[Log[x+I y]],{x,-4,4},{y,-4,4}] (*visualization of branch cut of complex function*)

plotcomplex[f_,N_]:=Module[{fn=f/.z->x+I y},Show[ContourPlot[Im[fn],{x,-N,N},{y,-N,N},ContourShading->Automatic,ExclusionsStyle->Red],ContourPlot[Re[fn],{x,-N,N},{y,-N,N},ContourShading->False,ContourStyle->Blue,ExclusionsStyle->Red],FrameLabel->{"Re(z)","Im(z)"},Background->Lighter[Orange],PlotRangePadding->0,PlotLabel->Framed[Grid[{{Style["-",Bold,Blue],"Real part"},{Style["-",Bold,Gray],"Imaginary part"}},Alignment->Left],FrameStyle->None,Background->White,RoundingRadius->5],ImageSize->300]] (*2D plot in complex plane with branch cuts in red*)

GraphicsRow[{plotcomplex[Log[z],10],plotcomplex[Log[z,π/4],10],plotcomplex[z*Tanh[π z]*Log[z^2+1],10]}] (*textbook examples*)

complexGrid=Compile[{{max,_Real},{n,_Integer}},Block[{r},r=Range[-max,max,2 max/(n-1)];
Outer[Plus,-I r,r]]];

complexHSB=Compile[{{Z,_Complex,2}},Block[{h,s,b,b2},h=Arg[Z]/(2 Pi);
s=Abs[Sin[2 Pi Abs[Z]]];
b=Sqrt[Sqrt[Abs[Sin[2 Pi Im[Z]] Sin[2 Pi Re[Z]]]]];
b2=0.5 ((1-s)+b+Sqrt[(1-s-b)^2+0.01]);
Transpose[{h,Sqrt[s],b2},{3,1,2}]]];

domainImage[func_,max_,n_]:=ImageResize[ColorConvert[Image[complexHSB@func@complexGrid[max,2 n],ColorSpace->"HSB"],"RGB"],n,Resampling->"Gaussian"];

domainPlot[func_: Identity,max_: Pi,n_: 500]:=Graphics[{},Frame->True,PlotRange->max,RotateLabel->False,FrameLabel->{"Re[z]","Im[z]","Domain Colouring of "<>ToString@StandardForm@func@"z"},BaseStyle->{FontFamily->"Calibri",12},Prolog->Inset[domainImage[func,max,n],{0,0},{Center,Center},2` max]];

domainPlot[] (*plot untransformed complex plane to understand what the colours indicate: brightness indicates absolute magnitude, saturation represents Im and Re parts*)

domainPlot[Log]

domainPlot[(#^2+1)/(#^2-1)&] (*zeros are at ±I poles are at ±1*)

Compare with built in Mathematica function:  

A) ComplexPlot uses a cyclic color function over Arg[f] to identify features such as zeros, poles and essential singularities. The color function goes from  to  counterclockwise around zeros, clockwise around poles and infinite cycles near essential singularities

ComplexPlot[(z^2+1)/(z^2-1),{z,-2-2 I,2+2 I}] (*zeros are at ±I poles are at ±1*)

ComplexPlot[(z^2+1)/(z^2-1),{z,-2-2 I,2+2 I},PlotLegends->Automatic] (*include a legend showing how the colors vary from -π to -π*)

ComplexPlot[(z^2+1)/(z^2-1),{z,-2-2 I,2+2 I},ColorFunction->"CyclicLogAbsArg"] (*use color shading to highlight features of the function*)

B) ComplexPlot3D plots Abs[f] with a cyclic color function over Arg[f] to identify features such as zeros, poles and essential singularities. The color function goes from  to  counterclockwise around zeros, clockwise around poles and infinite cycles near essential singularities:

ComplexPlot3D[(z^2 + 1)/(z^2 - 1), {z, -2 - 2 I, 2 + 2 I}](*zeros are at ±I poles are at ±1*)

ComplexPlot3D[(z^2+1)/(z^2-1),{z,-2-2 I,2+2 I},PlotLegends->Automatic](*include a legend showing how the colors vary from -π to -π*)

ComplexPlot3D[(z^2+1)/(z^2-1),{z,-2-2 I,2+2 I},ColorFunction->"CyclicLogAbsArg"](*use color shading to highlight features of the function*)

7. Calculus

7.1 Derivatives

D[x Sin[x],x]

(∂_x)(x Cos[x]+Sin[x])

fermi=1/(Exp[(e-ef)/(kbT)]+1) (*Fermi-Dirac distribution function from quantum statistical mechanics of noninteracting fermions*)

domainPlot[fermi/.{e->#,ef->1,kbT->0.2}&](*Poles of Fermi function in the complex plain are at the so-called Matsubara frequencies z=i*(2n+1)*Pi*kB*T *)

dfermi=-D[fermi,e] //FullSimplify(*Derivative of Fermi-Dirac function; add //FullSimplify to derive expression suitable for coding*)

Manipulate[Plot[{fermi,dfermi}/.{kbT->k,ef->fermienergy},{e,0,10}],{k,0.01,1},{fermienergy,0.01,10}]

bose=1/(Exp[(e-mu)/(kbT)]-1)(*Bose-Einstein distribution function from quantum statistical mechanics of noninteracting bosons*)

Manipulate[Plot[bose/.{kbT->k},{e,0,10}],{k,1,10},{mu,-10,0}]

D[x^3 z+2 y^2 x+y z^3,y,z] (*multivariable derivative*)

7.2 Integrals

Integrate[p x/(q+r x),x]

D[%,x]

% // Simplify

Integrate[x^2* Exp[x],{x,0,1}]

Integrate[Sin[Sin[x]],{x,0,1}]

NIntegrate[Sin[Sin[x]],{x,0,1}]

NIntegrate[Exp[-10^8 x^2],{x,-1,1}]

Plot[Exp[-10^8 x^2],{x,-0.1,0.1}]

NIntegrate[Exp[-10^8 x^2],{x,-1,1},WorkingPrecision->32,MaxRecursion->20] (*Use more adaptive recursion to resolve the integral of a rapidly varying function*)

Integrate[x* Sin[y^2]+y *Cos[x^2],{x,-1,1},{y,-2,x}] (*multivariable integral*)

N[%,10]

7.3 Infinite series

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

zex=Function[{g},Integrate[Exp[-x^2-g*x^4],{x,-Infinity,Infinity},Assumptions->g≥0]]

zser=Function[{g,n},Sum[(-g)^k/k!*Integrate[x^(4*k)*Exp[-x^2],{x,-Infinity,Infinity}],{k,0,n}]]

zex[1/50]

N[%,20]

N[zser[1/50,n=10],20]

Timing[Plot[zser[g=1/50,n],{n,0,50}]]

Timing[Plot[zser[g=1,n],{n,0,50}]] (*but optimal truncation stops working for "large" coupling constant g*)

7.4 Taylor series and Pade approximants

yam=Log[1+z]/z (*has branch point at z=-1*)

yamtaylor=Normal[Series[yam,{z,0,4}]] (*truncated Taylor series, to order N=4 around z=0, converges only within the circle |z|<1 and deviates from the exact function yam for|z|>1*)

yampd=PadeApproximant[yamtaylor,{z,0,{2,2}}] (*Taylor series is re-arranged into Pade [2/2] approximant*)

domainPlot[((1+((7 #)/10)+(((#)^2)/30))/(1+((6 #)/5)+((3 ((#)^2))/10)))&]

Plot[{yam,yamtaylor,yampd},{z,0,15},PlotRange->{{0,2}},PlotLegends->"Expressions"]

LESSONS LEARNED: Pade [2/2] well-approximates the original function yam with very high precession even beyond the radius of convergence up to Re z ∼ 10. Therefore, divergent power series  (Taylor series)  still contains information about the original function outside the convergence radius, and rearranging the coefficients of the expansion into the Pade approximants recovers this information. As a result, the conversion from the Taylor-form to the Pade-form usually accelerates the convergence and often allows good accuracy even outside the radius of convergence of the power series.

7.5 Differential equations

DSolve[y'[x]+y[x]==1,y[x],x]

y[x]+2*y'[x]+y[0]/.% (*solution as a rule applies only to y[x] and not to y[0] or y'[x]*)

DSolve[y'[x]+y[x]==1,y,x] (*gives solution as pure function*)

y[x]+2*y'[x]+y[0]/.%  (*now the rule applies to all occurences of y*)

y'[x]+y[x]==1/.%% (*substite solution into the original equation*)

DSolve[{x'[s]==Cos[t[s]],y'[s]==Sin[t[s]],t'[s]==s,x[0]==0,y[0]==0,t[0]==0},{x,y,t},s] (*Differential equation with initial conditions*)

ParametricPlot[Evaluate[{x[s],y[s]}/. %],{s,-10,10}] (*Parametric plot of the solution is Cornu spiral*)

7.6 Green’s functions (solving differential equations with delta function source)

circ=L i'[t]+R i[t] (*differential operator for the circuit where a voltage source v[t] is connected to a resistor R and an inductor L*)

gf=GreenFunction[{circ,i[0]==0},i[t],{t,0,∞},s] (*Green's function for this operator*)

v[t_]:=Sin[ω t]; Integrate[v[s] gf,{s,0,t},Assumptions->t>0] (*find current i[t] for a given voltage source v[t]*)

Plot[%/. {L->2,R->1,ω->5},{t,0,10}] (*visualize solution from previous cell*)

8. Special Functions of Mathematical Physics

Gamma                           	Gamma[z], Gamma[a,z]
Error                                 	Erf[z], Erfc[a,z], Erfi[z]
Fresnel integrals              	FresnelS, FresnelC
sin, cos integrals              	SinIntegral, CosIntegral
Legendre                         	LegendreP[n,x], LegendreQ
spherical harmonics         	SphericalHarmonicY[λ,m,θ,ϕ] 
Bessel                              	BesselJ[n,z],BesselY
modified Bessel          	BesselK[n,z],BesselI
Laguerre                          	LaguerreL
Hermite                              HermiteH[n,x]
hypergeometric               	Hypergeometric2F1, HypergeometricPFQ

Plot[BesselJ[2,x], {x,0,20}]

Manipulate[HermiteH[n,x],{n,0,20,1}]

si=Integrate[Sin[Exp[x]],x]

D[si,x]

si/.x->1.

Plot[si,{x,0,3}]

A=Cos[2*Pi*t/T]

J=Function[n,1/T*Integrate[Exp[I*(A-n*2*Pi*t/T)],{t,0,T}]] (*Integral appearing in the theory of Floquet topological insulators*)

J[5]

?BesselJ

9. Linear Algebra

Transpose[m]	                     transpose mᵀ
ConjugateTranspose[m]	conjugate transpose m† (Hermitian conjugate)
Inverse[m] 	                      	matrix inverse
Det[m]	                               		determinant
Minors[m]	                     	matrix of minors
Minors[m,k]	                    		k^(th) minors
Tr[m]	                               		trace
MatrixRank[m]	          	rank of matrix

$Post=If[MatrixQ[#],MatrixForm[#],#]& (*show matrices in human readable form*)

sigmax={{0,1},{1,0}}

sigmay={{0,-I},{I,0}}

sigmaz={{1,0},{0,-1}}

Table[PauliMatrix[k]//MatrixForm,{k,0,3}]

Heisenberg=KroneckerProduct[sigmax,IdentityMatrix[2]].KroneckerProduct[IdentityMatrix[2],sigmax]+KroneckerProduct[sigmay,IdentityMatrix[2]].KroneckerProduct[IdentityMatrix[2],sigmay]+KroneckerProduct[sigmaz,IdentityMatrix[2]].KroneckerProduct[IdentityMatrix[2],sigmaz] (*Hamiltonian of exchange interaction between two spins 1/2 or, equivalently, qubits*)

Tr[Heisenberg]

Heisenberg==ConjugateTranspose[Heisenberg]

A={{a,b},{c,d}}

B={{e},{f}}

KroneckerProduct[A,B] (*how Kronecker product of matrices works*)

KroneckerProduct[B,A] (*commuting order of matrices gives different Kronecker product*)

KroneckerProduct[sigmax,IdentityMatrix[2]].KroneckerProduct[IdentityMatrix[2],sigmax]== KroneckerProduct[sigmax,sigmax](*check theorem (A kron B).(C kron D)=(A.C)kron(B.D)*)

Eigenvalues[Heisenberg]

Eigenvectors[Heisenberg]

Heisenberg.{0,-1,1,0}==-3*{0,-1,1,0}