Mark G. Alford, Washington University Physics Department.
Revised 2020-08-24
This is an introduction to the most basic usage of Mathematica.
For more advanced examples, see Mathematica Techniques.
Full documentation is available from The Mathematica website including a Fast Introduction. There are tutorials at Mathematica Resources.
mathYou may need to type <<JavaGraphics.m to initialize plotting. The command-line interface allows convenient input editing: you can scroll back through previous commands using the up and down arrows, and edit them for re-use. For details see the text-based interface documentation.
In Mathematica, each command you type yields a result. If you terminate the command with a semicolon then the result will discarded and not printed out. This allows you to put several commands on the same line. Much of the notation will be familiar if you know how to program in C or Python.
Comments are enclosed in parentheses with asterisks (* ... *). Comments can be nested.
You can write mathematical expressions using the usual operators:x = 2*(3+4/2) 10 x^2 100 y = 2 (3+4/2) (* you can leave out the multiplication operator *) 10 1 / 2Pi (* but be careful to use parentheses where necessary *) Pi / 2 1 / (2 Pi) 1/(2*Pi) x y (* same as x*y *) 100 x=13; y=17; x*y (* multiple expressions on one line. Only the value of the last expression is "delivered" as the result *) 221Mathematica treats numbers without decimal points as exact, and does not numerically evaluate them. To force numerical evaluation, use the N function, or put a decimal point in a number somewhere:
x = 70/10 7 x = 70/9 (* No decimal point, so not evaluated *) 70/9 x = 70.0/9 7.77778 x = N[70/9] 7.77778Exponential notation for numbers (e.g. "3.5e15" in C) is done by writing "*^",
x = 3.5*^15;
x = 3.5 * 10^15; (* same thing *)
To undo assignments and make x and y revert to
being symbols of unknown value,
Clear[x]; Clear[y];
Sqrt[2]*Cos[Pi]*Exp[-1] -Sqrt[2]/E ArcTan[4] ArcTan[4] (* No decimal points, so not evaluated! *) ArcTan[4.0] 1.32582 (* all trig functions are in radians *) 5! (* factorial *) 120 Gamma[6] (* gamma function *) 120
phi[x_] = 1/(1+x^2); (* Now try using this function in various ways *) phi[0.0] 1.0 phi[1.5] 0.307692 phi[10.0] 0.00990099 phi[z] (1+z^2)^(-1) phi[ Sqrt[w] ] (1+w)^(-1)We can define a generalized version genphi whose width w must also be specified,
genphi[x_,w_] = 1/(1+(x/w)^2); genphi[0.0,0.2] 1.0 genphi[0.8,0.2] (* now x>>w so the function has a small value *) 0.0588235 genphi[x,5] (1 + x^2/25)^(-1)For serious programming, you will want to create functions that have local variables. This is done using the Module construction. Here is an example, with local variables sx and sy:
bump$fn[x_,y_,a_] = Module[{sx,sy}, (* declare local variables sx,sy *) sx = x/a; sy = y/a; 1/(1+sx^2 + sy^2) (* returned value: don't put a semicolon after it or you'll get "Null"! *) ];You can also exit from the function and return a result res by using Return[res];.
In these examples we defined the function using the = sign
which is the immediate evaluation
symbol, specifying that
the right hand side is to be evaluated immediately.
We will see in Mathematica Techniques
that for functions that involve numerical calculations
one must use the delayed evaluation
symbol := which evaluates the right hand side only when
the function is called.
Simplify[ Cos[theta] Cos[phi] - Sin[theta] Sin[phi] ] (* no assumptions needed *) Cos[phi + theta] Simplify[ Sqrt[x^2] * ArcCos[Cos[y]], {x>0,y>0,y<Pi} ] (* assumptions needed *) x yFor expressions involving logarithms or trig functions, use FullSimplify which has a more powerful algorithm. It may still need to be told what assumptions to make:
Simplify[ ArcSinh[x] - Log[x + Sqrt[1+x^2]], {x>0} ] ArcSinh[x] - Log[ x + Sqrt[1+x^2] ] (* no luck *) FullSimplify[ ArcSinh[x] - Log[x + Sqrt[1+x^2]] ] ArcSinh[x] - Log[ x + Sqrt[1+x^2] ] (* no luck *) FullSimplify[ ArcSinh[x] - Log[x + Sqrt[1+x^2]], {x>0} ] 0 (* at last! *)
Plot[ Cos[x], {x,-Pi,Pi} ];
plot1 = Plot[ Cos[x], {x,-Pi,Pi} ];Showing multiple functions on the same plot:
Plot[ {Sin[x],Cos[x]}, {x,-Pi,Pi} ];
Plot[ Cos[x], {x,-Pi,Pi}, PlotRange->{-0.5,0.5}];To control the positioning of the axes, use the AxesOrigin option:
Plot[ Sqrt[1+x^2], {x,-3,3} ]; (* y-axis begins at y=1 *) Plot[ Sqrt[1+x^2], {x,-3,3}, AxesOrigin->{0,0} ]; (* Axes cross at (0,0) *)To make the plot window larger on your screen, use the ImageSize option:
Plot[ Cos[x], {x,-Pi,Pi}, ImageSize->Scaled[2] ];To write the plot to PDF file, create it as a Graphics object and then "Export" it:
plot1 = Plot[ Cos[x], {x,-Pi,Pi} ]; Export["file.pdf",plot1];Mathematica infers the format from the dot-extension of the filename. Other possibilities include png, eps, svg, etc.
pl3d = Plot3D[ 1/(x^2+y+1), {x,0,10}, {y,0,20} ];
Integrate[ 1-x^2 ,{x,0,1}]
2/3
The limits can be symbols:
Integrate[ Cos[x], {x,a,b}]
-Sin[a] + Sin[b]
For a higher-dimensional integral, just give more integration variables
along with their limits:
Integrate[ Cos[x+y], {x,0,1}, {y,-a,a}]
2 Sin[1] Sin[a]
Sometimes, in order to get the desired answer you
need to give information about the parameters in the form of an
"Assumptions" option:
Integrate[ r^2*Exp[-r/a],{r,0,Infinity}, Assumptions->{a>0} ]
2 a^3
NIntegrate[ 1-x^2, {x,0,1}]
0.666667
or, integrate a user-defined function,
phi[x_] = 1/(1+x^2); NIntegrate[ phi[x],{x,0,Infinity}] (* mathematica can integrate to infinity *) 1.5708 genphi[x_,w_] = 1/(1+(x/w)^2); NIntegrate[ genphi[x,y],{x,0,Pi},{y,-1,1}] (* two-dimensional numerical integral *) 1.36271
Exp[I*3.0]
-0.989992 + 0.14112 I
You can obtain real and imaginary parts of any expression using
the Re and Im functions.
Functions can accept complex arguments, so you can write things like
Cos[I*x] Cosh[x] NIntegrate[ Exp[3*I*x] * Exp[-(x-1)^2],{x,-Infinity,Infinity}] -0.184946 + 0.0263634 IComplex conjugation is done by the Conjugate function:
Conjugate[ 3 + 4 I ]
3 - 4 I
You cannot plot a complex function, but
you can plot the real and/or imaginary parts of a complex function:
f[x_] = Exp[(-2+3*I)*x]; Plot[ Re[f[x]], {x,0,2}]; Plot[ {Re[f[x]],Im[f[x]]}, {x,0,2}];
v = {2.0, -1.0, 0.0} ;A matrix or two-dimensional array is a list of lists:
M = { {1.3, -2.0, 0.9}, {-2.2, 0.3, 1.1}, {-1.4, 1.5, -0.2} };Matrix multiplication or dot product uses the dot operator. Scalar multiplication uses "*",
w = 0.5 * M.v
{2.3, -2.35, -2.15}
w.v
6.95
There are many useful built-in matrix functions:
Minv = Inverse[M]; MT = Transpose[M]; evals = Eigenvalues[M] {2.91742, -1.69075, 0.173336} evecs = Eigenvectors[M] {{-0.512376, 0.661041, 0.548175}, {-0.587781, -0.775991, 0.228806}, {0.397766, 0.553388, 0.731809}} Tr[M] (* Note Trace is "Tr", not "Trace" *) 1.4You can even take the exponential of a matrix,
Mexp = MatrixExp[M];Of course, you can calculate the norm of a vector or matrix,
Norm[v] (* sqrt of sum of absolute value squared of elements *) 2.23607 Norm[M] 3.59095User-defined functions of matrices are just like any other function. Here are some useful ones:
Commutator[a_,b_] = a.b - b.a; Anticommutator[a_,b_]= a.b + b.a; Adjoint[m_] = Transpose[Conjugate[m]];Subscripting uses double square brackets. Subscripts start at 1. This is one of the differences between Mathematica and C or Python. (The "zeroth" element of a list is its type, actually).
v[[3]] 0. M[[2,3]] (* treat M like a matrix *) 1.1 M[[2]][[3]] (* treat M like a list of lists *) 1.1Printing matrices in a more readable form:
Print[MatrixForm[M]];
f[x_]=Exp[-x^2]; f'[x] (* first derivative, prime notation *) -2*x*E^(-x^2) f''[x] (* second derivative, prime notation *) -2*E^(-x^2) + 4*x^2*E^(-x^2) D[f[x],{x,5}]; (* fifth derivative, using D function *) Derivative[5][f][x0]; (* fifth derivative of f, at x=x0 *)Function of more than one variable:
f[x_,y_] = Exp[-x^2+x*y]; D[ f[x,y], x] (* derivative with respect to x, at fixed y *) E^(-x^2 + x*y)*(-2*x + y) D[ f[x,y], {y,2} ] (* second derivative wrt y, at fixed x *) E^(-x^2 + x*y)*x^2 Derivative[5,2][f][x0,y0]; (* fifth deriv wrt x and second deriv wrt y, evaluated at (x0,y0) *)
Series[ Cos[x], {x,0,4}]
1 - x^2/2 + x^4/24 + O[x]^5
Use the Normal function to strip off the + O[x]^5
and get a regular polynomial:
Normal[Series[ Cos[x], {x,0,4}]]
1 - x^2/2 + x^4/24
Expand Log[x] around x=1 to 3rd order:
Normal[Series[ Log[x], {x,1,3}]]
-1 + x - (-1 + x)^2/2 + (-1 + x)^3/3
The Series function will often produce a simpler answer when you
give it some helpful information about the parameters by using the "Assumptions" option:
Normal[Series[ Sqrt[ x^2 - a^2 ] , {x,a,1}, Assumptions->{a>0} ]]
Sqrt[2] Sqrt[a] Sqrt[x-a]
Without the assumption that a is real and positive,
Mathematica would deliver a much more complicated-looking result.
Normal[Series[ Log[1+epsilon], {epsilon,0,3}]]
epsilon - epsilon^2/2 + epsilon^3/3
Print["The log of ",3.5*^15," is ",Log[3.5*^15] ]; 15 The log of 3.5 10 is 35.7915 Print["The log of ",3.5*^15," is ",SetPrecision[Log[3.5*^15],8] ]; 15 The log of 3.5 10 is 35.791539If you don't like the way exponentials are printed, use the "CForm" function to get the more standard "e" notation:
Print[CForm[3.5*^15]]; 3.5e15 Print[CForm[ SetPrecision[Pi*10^10,5] ]]; 3.1416e10
For more advanced examples, see Mathematica Techniques.