Mark G. Alford, Washington University Physics Department.
Revised 20170324
For an introduction to the basics of Mathematica, see my
Introduction to Mathematica.
Full documentation is available from The Mathematica website.
For serious projects you will need to create "package" files containing definitions of the functions that perform your calculations. Here are two ways to do it.
<< relaxation_solver.mMathematica searches for this file in the directories (folders) listed in the variable $Path. To make sure Mathematica finds your file, you'll need to include in $Path the directory where the file is located. E.g., if it is in a folder math that is a subdirectory in your home directory,
AppendTo[$Path, FileNameJoin[{$HomeDirectory, "math"}]] << relaxation_solver.m
M$of$I$sphere[radius_,density_]=Module[{r$axial,z}, (* declare r$axial and z as local symbols, internal to the function *) r$axial[z_]=Sqrt[radius^2z^2]; (* you can define a function within a function *) Integrate[ Pi/2 * density * r$axial[z]^4, {z,radius,radius}] ]; M$of$I$sphere[R,rho] (* use the function *) (8*Pi*rho*R^5)/15All the symbols that we use temporarily inside the module (r$axial and z in this case) are declared at the start of the module as local variables. This prevents interference with any other entity called z or r$axial. We do not declare radius or density as local because they are arguments to the function we are defining.
"a and b"  a && b  
"a or b"  a  b  
"not a"  !a 
First example:
testing whether a list contains a given element.
We define
the function using delayed evaluation (:=
) because it cannot
be evaluated until it is called with an actual list and element as arguments.
contains[list_,element_]:=Module[{i,found}, (* all variables internal to the function are declared as local *) i=1; (* Loop counter *) found=False; (* 'found' is a Boolean variable *) While[i<=Length[list] && ! found, (* Repeat the loop as long as the search has not succeeded and we haven't reached the end of the list *) found= (list[[i]]==element); (* Check for a match: As in C and Python, we useSecond example: counting how many times a given element appears in a list:==for testing the equality of two quantities *) i+=1; (* Increment the loop counter: As in C and Python, i+=1 means i=i+1*) ]; found (* returns a Boolean *) ]; contains[{4,9,2,8,1,3}, 3] (* Use the function we have just defined *) True
occurrences[list_,element_]:=Module[{i,nfound}, nfound=0; For[i=1, i<=Length[list], i++, (* As in C and Python, i++ means i=i+1 *) If[ list[[i]]==element, nfound+=1]; ]; If[nfound>0, Print["found ",nfound," occurrences"] ,(*else*) Print["found no occurrences"] ]; nfound (* returns number of occurrences *) ];
eqn = ( x^2  c == 0 ); (* As in C and Python, we useWe see here that the Solve command returns a list of solutions. Each solution is itself a list of substitution rules, taking the form x > EXPRESSION. To set some variable x$soln equal to the solution you have to apply the relevant set of substitution rules to the unknown variable that was in the equation (x in this case).==for testing the equality of two quantities *) rules$list = Solve[ eqn, x] {{x > Sqrt[c]}, {x > Sqrt[c]}} (* There are two solutions in the list *) x$soln = x /. rules$list[[2]] (* to get the second solution, apply the second set of rules to x *) Sqrt[c])
twovar$root[c_] = Module[{x,y,eqn$list,rules$list}, (* all variables internal to the function are declared as local *) eqn$list = { x^2 +y^2 == 2*c, x + y == 0}; rules$list = Solve[ eqn$list, {x,y} ]; (* solve the two equations for x and y *) (* At this point, rules$list is {{x > Sqrt[c], y > Sqrt[c]}, {x > Sqrt[c], y > Sqrt[c]}} *) {x,y} /. rules$list[[2]] ]; twovar$root[w] (* use the function we just defined *) {Sqrt[w], Sqrt[w]}
list1 = {1,2,3,4,5,6,7,8,9};
list1[[4;;8]]
{4, 5, 6, 7, 8}
The index 1 means the last element of the list, 2 means the nexttolast,
and so on:
list1[[5;;1]] {5, 6, 7, 8, 9} list1[[4;;1]] {6, 7, 8, 9}Lists can be built using the Table function:
list2 = Table[ i, {i,1,9} ]
{1, 2, 3, 4, 5, 6, 7, 8, 9}
To populate lists with more complicated contents, we can use the fact that
anywhere a value is required, one can put a series of statements,
separated by semicolons, that yield
that value. For example, to make a series of {x,f[x]} pairs:
xy$list = Table[x=i/5.0; y=Tanh[x]; {x,y} , {i,0,10} ]
{{0., 0.}, {0.2, 0.197375}, {0.4, 0.379949}, {0.6, 0.53705},
{0.8, 0.664037}, {1., 0.761594}, {1.2, 0.833655}, {1.4, 0.885352},
{1.6, 0.921669}, {1.8, 0.946806}, {2., 0.964028}
}
LogPlot[ 1.5^x , {x,0,10} ] (* logarithmic y axis *) LogLinearPlot[ Tanh[x] , {x,0.01,10} ] (* logarithmic x axis *) LogLogPlot[ x^1.67, {x,0.01,10} ] (* logarithmic x and y axis *)One can plot lists of (x,y) coordinates instead of functions:
xylist = Table[ x=i/10.0; y=x^2; {x,y}, {i,1,10} ] (* 10 points on the y=x^2 curve *) {{0.1, 0.01}, {0.2, 0.04}, {0.3, 0.09}, {0.4, 0.16}, {0.5, 0.25}, {0.6, 0.36}, {0.7, 0.49}, {0.8, 0.64}, {0.9, 0.81}, {1., 1.}} ListPlot[xylist];Lists can also be plotted on logarithmic axes:
ListLogPlot[xylist]; (* y axis is a log axis *) ListLogLogPlot[xylist]; (* both axes are log *)To plot functions and lists on the same plot, use the Show function to combine the plots. In this case, when we create plot1 and plot2 we stop them from being displayed by setting the DisplayFunction to do nothing (
Identity), and then set it back to the default ($DisplayFunction) when we combine the plots using Show:
plot1 = Plot[ Cos[x*Pi/2], {x,1,1}, DisplayFunction>Identity ]; xylist = Table[ x=i/10; y=1x^2; {x,y}, {i,1,10} ]; plot2 = ListPlot[xylist,DisplayFunction>Identity ]; plot$12 = Show[plot1,plot2]; (* combine plot1 and plot2 on a single plot *) Show[plot$12, DisplayFunction>$DisplayFunction]; (* show the combined plot on the screen *) Export["plots_comb.pdf",plot$12] (* create PDF file of the combined plots *)To control color, thickness, etc, use the PlotStyle option:
(* one plot style for one curve *)
Plot[ Cos[x*Pi/2], {x,1,1},
PlotStyle>{Cyan,Thickness[0.007],Dashing[{0.03,0.03}]}
];

(* list of plot styles for different curves *)
Plot[ {Cos[x*Pi/2],1x^2*Pi^2/8}, {x,1,1}, PlotStyle>{
{Cyan,Thickness[0.01],Dashing[{0.03,0.05}]},
{Purple, Thickness[0.005]}
} ];

x$rule = FindRoot[ x^3 Tanh[x], {x,0,1.0} ] {x > 0.} x$soln = x /. x$rule 0.As with Solve, the result is returned in the form of a rule, but FindRoot returns at most one solution, which in this case may not be what was intended: there are two roots in the range x=0 to 1, and FindRoot found the trivial one. Specifying a narrower initial search range helps it find the nontrivial solution:
x$rule = FindRoot[ x^3 Tanh[x], {x,1.0,1.2} ];
x$soln = x /. x$rule
0.893395
Note that Mathematica will go outside the initial search range
(in this case x=1.0 to 1.2) if necessary.xy$solutions[x$guess_?NumericQ, y$guess_?NumericQ] := Module[{x,y,xy$rule}, (* we will explain below why the?NumericQpatterns are needed *) xy$rule = FindRoot[ {x^3Tanh[y], y  x^2}, { {x,x$guess}, {y,y$guess} } ]; {x,y} /. xy$rule ]; xy$solutions[0.7,0.5] {0.853834, 0.729033}
eqn = ( D[theta[t],{t,2}] == theta[t] ); (* or we could have written it as theta''[t] == theta[t] *) bc = {theta[0]==1, theta'[0]==0}; (* boundary conditions *) soln$rules$list = DSolve[ {eqn, bc} , theta[t], t] {{theta[t] > Cos[t]}} theta$soln[t_] = theta[t] /. soln$rules$list[[1]] (* there is only one solution, but we still have to specify the first solution *) Cos[t]DSolve returns a list (of length 1) containing a list of substitution rules (see Rules... above), one for each unknown function in the differential equation. In this case there is one unknown function theta, so soln$rules$list[[1]] is a list containing one rule that replaces instances of the unknown theta[t] with the solution to the differential equation.
theta$soln''[t] + theta$soln[t]
0
SHO$solution[tmax_?NumericQ]:= Module[{eqn,bc,theta,t,soln$rules$list}, (* we will explain below why theLike DSolve, NDSolve returns a list (of length 1) containing a list of substitution rules, one for each unknown function in the differential equation. In this case soln$rules$list[[1]] is a list containing one rule that replaces instances of the unknown theta with the numerical solution to the differential equation in the form of an?NumericQpattern is needed *) eqn = ( D[theta[t],{t,2}] == 2*theta[t] ); (* or theta''[t] == 2*theta[t] *) bc = {theta[0]==0.5, theta'[0]==0}; (* boundary conditions *) soln$rules$list = NDSolve[ {eqn, bc} , theta, {t,0,tmax}]; (* We evolve from t=0 to tmax *) (* At thus point soln$rules$list is {{theta > InterpolatingFunction[{{0., }}, <>]}} *) theta /. soln$rules$list[[1]] (* NDSolve always gives one solution, but you still have to ask for the first solution *) ]; phi$IF = SHO$solution[3.14]; (* Try using the function we have just defined *) (* phi$IF is now an InterpolatingFunction, which can be evaluated, plotted, etc *) phi$IF[0] 0.5 phi$IF[2.5] 0.461702 Plot[ phi$IF[tau], {tau,0,3.14}];
InterpolatingFunction. To get that Interpolating Function we have to apply that rule to theta. Mathematica treats the Interpolating Function like a numerically defined function, so you can plot it, numerically integrate it, etc.
fn$IF = Interpolation[ Table[ x=i/10.0; y=Sqrt[x]; {x,y}, {i,0,20}], InterpolationOrder>2]; (* Construct our own InterpolatingFunction *) Plot[{ fn$IF[x],Sqrt[x]},{x,0,2}]; (* interpolation works well except near x=0 *) (* Import the module that can look inside the InterpolatingFunction: *) Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]; (* We can obtain the x limits over which the interpolating function is defined, allowing us to plot it even if we didn't generate it ourselves *) {x$min,x$max} = InterpolatingFunctionDomain[fn$IF][[1]] {0., 2.} Plot[ fn$IF[x], {x,x$min,x$max} ]; (* We can obtain the xy points between which the function is interpolated *) xvals = InterpolatingFunctionCoordinates[fn$IF][[1]] {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.} yvals = InterpolatingFunctionValuesOnGrid[fn$IF] {0., 0.316228, 0.447214, 0.547723, 0.632456, 0.707107, 0.774597, 0.83666, 0.894427, 0.948683, 1., 1.04881, 1.09545, 1.14018, 1.18322, 1.22474, 1.26491, 1.30384, 1.34164, 1.3784, 1.41421} (* we can find out what order polynomial is used in the interpolation *) order = InterpolatingFunctionInterpolationOrder[fn$IF][[1]] 2 (* we can switch the x and y data points to make a new InterpolatingFunction that is approximately the inverse of the original function *) fn$inv$IF = Interpolation[ Transpose[ {yvals,xvals} ], InterpolationOrder>2 ]; fn$inv$IF[fn$IF[3.0]] 2.95755 (* The inversion is approximate *)
trans$root[a_] = Module[{x},
x /. Solve[ x^3  Tanh[x/a]==0, x]
];
Solve::nsmet: This system cannot be solved with the methods available to Solve.
OK, so we are forced to write it as a function that finds the answer numerically:
trans$root$numeric$naive[a_] := Module[{x}, x /. FindRoot[ x^3  Tanh[x/a]==0, {x,1} ] ];In this case we defined the function using delayed evaluation (
:=) because it can only be evaluated once a number has been provided as the value of the argument a. The problem is, this function doesn't "know" that its argument must be a number. Errors will occur if we later use this function in a context where Mathematica tries to feed a symbolic value to the function. For example,
NIntegrate[ trans$root$numeric$naive[a], {a,1/2,1} ]
FindRoot::nlnum: The function value {1.  1. Tanh[1./a]} is not a list of numbers with dimensions {1} ...
This surprising error arises because Mathematica, even though it
has been asked to
evaluate the integral numerically, attempts some symbolic manipulation
of the integrand, which is a purely numerical function.
To prevent this, we must define the function
so that it only accepts numbers as arguments.
To do this we add a _?NumericQ Pattern in the definition:
trans$root$numeric[a_?NumericQ] := Module[{x}, x /. FindRoot[ x^3  Tanh[x/a]==0, {x,1} ] ];NumericQ is a function that returns True if its argument is a number. a_?NumericQ is a pattern that says that a can only be replaced with a number. The function will now only accept numbers as arguments, and can be numerically integrated without error messages:
NIntegrate[ trans$root$numeric[a], {a,1/2,1} ]
0.472751
xydata = {{0., 0.951}, {0.1, 0.846}, {0.2, 0.309}, {0.3, 0.095}, {0.4, 0.694}, {0.5, 1.029}, {0.6, 0.706}, {0.7, 0.448}, {0.8, 0.768}, {0.9, 1.171}, {1., 0.963} };
nlmfit = NonlinearModelFit[ xydata, a + b*Cos[omega*t], {{a,0.2}, {b,1.0}, {omega,5}}, t]; (* we give an initial guess for each parameter *)
NonlinearModelFit can be rather sensitive to the initial guesses for the parameter values. If the fit is successful then
nlmfit is now a FittedModel object, which superficially behaves like the bestfit function to the data, but also contains a lot of information about the fit.
Normal[nlmfit] (* shows the best fit function *) 0.0707912 + 1.01908 Cos[6.44376 t] nlmfit[0.75] (* evaluate the best fit function at t=0.75 *) 0.193222 (* We can make a plot of the data and the fit: *) plot$data=ListPlot[xydata,DisplayFunction>Identity]; plot$fit=Plot[nlmfit[t],{t,0,1},DisplayFunction>Identity]; Show[{plot$fit,plot$data},PlotRange>{2,2},DisplayFunction>$DisplayFunction];
nlmfit["BestFitParameters"] {a > 0.0707912, b > 1.01908, omega > 6.44376} nlmfit["ParameterErrors"] {0.0621924, 0.0790336, 0.162831} (* a is not significantly different from 0, but b and omega are *) nlmfit["ParameterTable"] Estimate Standard Error tStatistic PValue a 0.0707912 0.0621924 1.13826 0.287949 b 1.01908 0.0790336 12.8943 1.23746*^6 omega 6.44376 0.162831 39.5732 1.82829*^10 nlmfit["AdjustedRSquared"] 0.940382We can obtain many other diagnostics and properties of the fit, including an ANOVA table. See the Mathematica documentation for full details.