Overview of Numerical Computations in Mathematica

Rob Knapp

Wolfram Research, Inc.

Introduction

The symbolic and numerical components of Mathematica are tightly integrated and interwoven in many ways.  The graph below gives a rough sketch of the ways different parts are interrelated.

ar = "Arithmetic";
ne = "Numerical Evaluation";
sc = "Symbolic Computation";
la = "Linear Algebra";
de = "Differential Equations";
op = "Optimization";
rt = "Rootfinding";
qd = "Quadrature";
LayeredGraphPlot[{ar->ne, ar->sc, ar->la, ar->de, ar->op, ar->rt, ar->qd, ne->sc, ne->de, ne->op, ne->rt, ne->qd, sc->de, sc->op, sc->rt, sc->qd, la->sc, la->de, la->op, la->rt, op->sc},DirectedEdges->True,  VertexLabeling->All]

In this presentation I plan to discuss each of these parts except for symbolic computation, starting from the top.

Arithmetic

Of course this sounds too simple to even bother with, but really, it is not.  Mark Sofroniou made a nice Wolfram Blog post on "Arithmetic Is Hard--To Get Right" recently.

Mathematica supports arithmetic for four different types of explicit numbers, Integer, Rational, Real, and Complex. Exact representations of numbers like (Sqrt[2])could be considered a fifth type, but those fall more under the purview of symbolic computation and numerical evaluation than arithmetic.

Internally, the four types are actually broken into subtypes -- the most apparent of these subtypes are machine reals and integers.  Mathematica uses machine hardware numbers when those are sufficient to represent the desired quantity to take advantage of the blinding speed of arithmetic on modern processors.

Integers too big or reals too precise to fit within the confines of machine numbers are still quite fast, in part thanks to our use of the Gnu Multiprecision Library (GMP).   Mark Sofroniou will be giving a presentation on "A Framework for Efficient Arbitrary-Precision Numerical Computation" later this afternoon.

In addition to arithmetic with like types, Mathematica supports a coercion model that allows you to evaluate expressions naturally without having to specify coercion.  This may seem like an obvious thing to you as Mathematica users, but not all systems support such a consistent model.  The key things behind the coercion model are

Exact in ⇒ exact out

The precision of the result is based on the precision of the inputs -- coercion is done to achieve the best precision justifiable from the inputs.

Precision of results is computed using significance arithmetic.

Machine numbers are considered to have less precision than any other numbers

When machine numbers are terms in an arithmetic operation, all the inputs will be coerced to machine numbers if possible.

This coercion model applies generally to computations in Mathematica, not just arithmetic.

Significance arithmetic keeps track of justified precision using a first order approximation.  It is not rigorous like interval arithmetic (also supported in Mathematica with Interval objects) but it is fast and, with reasonably high precision accurate as well.  The Precision command in Mathematica returns a measure of the relative size of the estimated error interval and the Accuracy command gives a measure of the absolute size of the error interval.

This gives the precision and accuracy for a number given with 30 decimal places.

x = 0.999999999999999999999999999999;

{Precision[x], Accuracy[x], Precision[32 x], Accuracy[32 x]}

For the number x above, the error interval is [ 0.9999999999999999999999999999985, 0.9999999999999999999999999999995], having length Δx  = (10^(-30)).

In this example, the result has no precision since all significant digits have cancelled.  The 0 in the output indicates that the error interval contains zero.

1 -x

Precision[%]

Most computations do not lose a lot of precision.  In fact, some even gain precision.

Repeatedly iterating the function 1 - x/2 gives a result with higher precision than the input.

Nest[1 - #/2&,0.999999999999999999999999999999, 1000]

For the example 1 - x/2, you can easily see why there is a precision gain:
Suppose that there is initially an error Δx (in this case about (10^(-30))).  Then

y + Δy = 1 - (x + Δx)/2 = 1 - x/2 - Δx/2

so the error in the result is about half that of the input.  Its slightly greater than that because of the arithmetic of subtracting from 1 (since in Mathematica numbers are internally represented in binary, division by 2 introduces no error).  Near 2/3, the error of the subtraction is less than the error Δx so you get about one bit extra per iteration.

x0 = 0.999999999999999999999999999999;
Show[ListPlot[Map[Precision, NestList[1 - #/2&, x0, 100]],PlotStyle->Red], Plot[Evaluate[Precision[x0] + k Log[10,2]],{k,1, 100}]]

The one other aspect of arithmetic I'll mention here is that it is threaded over lists in such a way that, for example adding two lists adds them elementwise.

This adds a list to twice another list.

{1,2,3} + 2 {-1,1,0}

This actually works with all the functions that have the Listable attribute:

listable=Select[Names["System`*"],MemberQ[Attributes[#],Listable]&]

In general, if you can use listability for operations on large lists, your computations will be much faster.

Compare adding a large nested random array to itself using listability to using indexing.

a = RandomReal[1,{100, 100, 100}];

Timing[aa1 = a + a;]

Timing[aa2 = Table[a[[i,j,k]]+a[[i,j,k]],{i,100},{j,100},{k,100}];]

aa1 === aa2

Numerical Evaluation

The command N in Mathematica  does numerical evaluation.  Generally when an expression has the NumericQ attribute it will evaluate to a number inside of N.  At a user level, this seems simple and natural, but there is a huge amount of machinery behind it.

Evaluate an integral to get a numeric expression representing the exact result.  Then find the numerical value with N.

in = Integrate[Sin[Exp[s]],{s,-1,1}]

{NumericQ[in], N[in]}

This example can be done because we have good numerical algorithms to compute the special sine integral function.  Mathematica supports a many special functions all of which require several algorithms for computing the results in different ranges of numbers.

The default for N is to use machine precision, using hardware numbers as much as possible.  However, N can work adaptively to any precision, so you can get the result of the integral as precise as you'd like.

This gives the result of the integral to 100 digits.

N[in,100]

Note that when functions, special or otherwise are evaluated, their conditioning is taken into account.  This way , you get a precision out justified by the arguments in.  The same thing happens for the hyperbolic function coth as for the example I showed before:

This iterates the hyperbolic function coth(x) 1000 times.

Nest[Coth, 1.11111111111111111111,1000]

Each time the function is evaluated, the precision of the result is what can be justified by the precision of the input.

One new way Mathematica has for representing numbers is via transcendental root objects:

trs = Reduce[Coth[x]==x,x,Reals]

These can be evaluated to any precision.

N[trs[[1,2]], 1000]

Linear Algebra

Mathematica supports linear algebra computations for machine numbers, higher precision numbers, as well as with exact arithmetic.

The machine number linear algebra is just as fast as in programs that are dedicated to linear algebra because it is based on a PackedArray representation of lists and uses a fast library for the Basic Linear Algebra Subroutines (BLAS) that is optimized on a processor basis.  The Intel mathkernel BLAS used on this machine is also able to use multithreading to access the dual core, giving additional speedup for big computations.

This computes the LU decomposition for a dense 1000x1000 matrix.

m = RandomReal[1,{2000, 2000}];
Timing[lu = LinearSolve[m]]

The decomposition is used in a functional form that can easily be used to solve linear systems.

b = m.Range[2000];
Timing[x = lu[b];]

Norm[(x - Range[2000])/Range[2000]]

For some matrices, machine number computations are inadequate.  For example, computing the smallest eigenvalue of a Hilbert matrix:

This assigns m to be the exact 100x100 Hilbert matrix (Subscript[a, i j])= 1/(i+j-1)

Short[m = HilbertMatrix[50]]

This computes a representation for the smallest eigenvalue.

Timing[ex = Eigenvalues[m, -1]]

This computes the numerical value using numerical evaluation.

N[ex]

This computes the value using machine number arithmetic.

Timing[em = Eigenvalues[N[m],-1]]

This computes the value using 75 digit arithmetic

Timing[eb = Eigenvalues[N[m, 75], -1]]

Note that the last example is done somewhat differently than numerical evaluation with N.  Instead it is done with a fixed precision based on the precision of the input matrix.  With the sheer number of operations required to do numerical linear algebra, using adaptive precision could lead to excess loss of precision.  Also, numerical linear algebra as a field is very well studied and understood with the assumption of fixed precision.

This shows that the result from using 75 digit arithmetic in linear algebra only has about 13 digits computed correctly.

(ex - eb)/ex

When it is more appropriate to use a fixed precision in a computation, you can do this by setting the system parameters $MinPrecision and $MaxPrecision to be the same.

This computes the positive root of coth(x)=x to 100 digits using fixed precision.

Block[{$MinPrecision = 100, $MaxPrecision = 100}, FixedPoint[Coth, 1`100]]

If you don't use fixed precision, the computation may terminate at a much higher precision than you need (or in some cases, may not terminate at all)

FixedPoint[Coth, 1`100]

FixedPoint[1-#/2&, 1`100, 1000]

Sparse Arrays

In the development process for Mathematica 5, some of the numerical functionality, in particular PDE solving and interior point methods for linear programming needed to use a sparse matrix representation along with sparse numerical linear algebra.

Though in the short term, it would have been easier to put together something specialized for these uses, we felt that, in the long term, the benefits of integrating them fully into the system was worth the extra work (and it was a lot of extra work)

SparseArrays also represent lists, vectors, and matrices, but the representation is made explicitly rather than hidden as in the Packed Array case.  Arithmetic and other operations typically work with SparseArrays, just as if it were a list, though.

Construct a tridiagonal matrix:

S = SparseArray[{{i_, i_}->-2, {i_, j_} /; Abs[i-j]==1->1.}, {10, 10}]

Show it in ordinary matrix form

MatrixForm[S]

Combine it with another matrix.

op = S/.1^2 + SparseArray[{{i_, i_}->1.},{10,10}]

Convert it to a normal matrix

Normal[op]

With SparseArray, linear algebra can work with extremely large matrices

Here is a relatively small sparse matrix:

S = ExampleData[{"Matrix","HB/bcsstm27"}]

ArrayPlot[S]

The LU decomposition is computed in a fraction of a second.

Timing[lu = LinearSolve[S];]

Numerical procedures working with functions.

The rest of what I'll cover today involves numerical procedures and algorithms that work with functions, such as numerically computing a definite integral for a function or finding the solution to a system of ordinary differential equations given a right hand side function.  Typically these methods work by sampling the function at points determined by the algorithm and error is estimated based on an approximation the procedure uses.

Though the different sorts of tasks typically done are quite different,  we have used several design principles to make the interfaces to these functions be as similar as possible.

Functions are task or problem based rather than algorithm based (e.g. FindRoot, NDSolve instead of Newton, VODE).

Particular algorithms and methods can usually be selected and/or tuned with Method options.

Otherwise symbolic and other analyses are used to try to determine an appropriate algorithm.

Output from one function can easily be used in other functions.

Specialized data structures (such as InterpolatingFunction returned from NDSolve) have interfaces that integrate
 them with the rest of the system.

Error is controlled by options AccuracyGoal and PrecisionGoal

Computations are typically done in fixed precision specified by WorkingPrecision.

For functions that iteratively improve the solution, a sequence of a fixed precisions may be used

EvaluationMonitor and StepMonitor options allow you to monitor the method being used.

Frameworks are set up that allow you to define and "plug-in" your own methods.

Differential Equations

KdV equation

Solving the KdV equation is an example that illustrates many of the design principles described above.

This sets up a display that will be used to monitor the progress of the solution.

Dynamic[plot]

This solves the KdV equation with NDSolve.

L = 10;
Timing[NDSolve[{D[u[t,x],t] == D[u[t,x],x,x,x] + 6 u[t,x]D[u[t,x],x],u[0,x]== Exp[-x^2],u[t,-L]==u[t,L]},u,{t,0,1},{x,-L,L},StepMonitor:>(plot = Plot[u[t,x],{x,-L,L}, PlotRange->{-1,1}])]]

Because of the large dispersion and with periodic boundary conditions, using a pseudospectral approximation for the derivative terms can give just as accurate a solution with far fewer spatial points.  This is controlled by a method option.

Solve the KdV equation with a pseudospectral method.

L = 10;
Timing[NDSolve[{D[u[t,x],t] == D[u[t,x],x,x,x] + 6 u[t,x]D[u[t,x],x],u[0,x]== Exp[-x^2],u[t,-L]==u[t,L]},u,{t,0,1},{x,-L,L},StepMonitor:>(plot = Plot[u[t,x],{x,-L,L}, PlotRange->{-1,1}]), Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->"Pseudospectral"}}]]

A lot of the time in the solution above is actually spent monitoring, since creating an intermediate solution and doing the plot is relatively expensive.

Use a discretized intermediate solution for the monitor.

L = 10;
Timing[NDSolve[{D[u[t,x],t] == D[u[t,x],x,x,x] + 6 u[t,x]D[u[t,x],x],u[0,x]== Exp[-x^2],u[t,-L]==u[t,L]},u,{t,0,1},{x,-L,L},StepMonitor:>(plot =ListPlot[u[t,x], PlotRange->{-1,1},DataRange->{-L,L}]), Method->{"MethodOfLines", "DiscretizedMonitorVariables"->True, "SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->"Pseudospectral"}}]]

This is fast enough that its now possible to enlarge the spatial domain.

Use a larger spatial domain.

L = 50;
Timing[NDSolve[{D[u[t,x],t] == D[u[t,x],x,x,x] + 6 u[t,x]D[u[t,x],x],u[0,x]== Exp[-x^2],u[t,-L]==u[t,L]},u,{t,0,1},{x,-L,L},StepMonitor:>(plot =ListPlot[u[t,x], PlotRange->{-1,1}, DataRange->{-L,L}]), Method->{"MethodOfLines", "DiscretizedMonitorVariables"->True, "SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->"Pseudospectral"}}]]

Ideally, for this computation, the best thing to do would be to solve up to the time at which the waves start to spill over the periodic boundary conditions.  Its quite easy to do this by considering that as an event.

Solve until the time at which the size of the solution at either boundary becomes above a threshold τ.

τ = 0.01;
Timing[sol = NDSolve[{D[u[t,x],t] == D[u[t,x],x,x,x] + 6 u[t,x]D[u[t,x],x],u[0,x]== Exp[-x^2],u[t,-L]==u[t,L]},u,{t,0,1000},{x,-L,L},StepMonitor:>(plot =ListPlot[u[t,x], PlotRange->{-1,1}, DataRange->{-L,L}]), Method->{"MethodOfLines", 
"DiscretizedMonitorVariables"->True, 
"SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->"Pseudospectral"},
Method->{"EventLocator", "Event":>(Max[Abs[u[t,x]][[{1,-1}]]] - 0.01), "EventAction":>Throw[end = t, "StopIntegration"]}}]]

And finally, the InterpolatingFunction returned is integrated with the system, so it is easy to make a nice visualization of the solution.

Plot3D[Evaluate[First[u[t,x] /. sol]],{t,0,end},{x,-L,L}, MaxRecursion->5, Mesh->False, PlotRange->All]

Numerical Integration

Numerical integration provides a good example of where symbolic preprocessing can help choose and assist the algorithms that are used.

Here is a discontinuous function.  An additional definition is given for a blackbox function with the same values, but effectively hidden definition.

f[x_] :=Sin[Floor[x]]/x
blackbox[x_?NumberQ] := f[x];
Plot[f[x],{x,0,10}, ExclusionsStyle->Directive[Red, Dashing[.005]]]

Compute the exact value of the integral for this example.

exact = Integrate[f[x],{x,0,10}]

Compute the value of the integral numerically.

nf = NIntegrate[f[x],{x,0,10}]

Check the error in the computed value.

nf - exact

NIntegrate has done extremely well computing the value numerically when it can se the definition of the function.

Compute the value of the same function with the definition hidden to NIntegrate, so that all it can do is sample values.

fb = NIntegrate[blackbox[x],{x,0,10}]

fb - exact

In the black box case, the result is not nearly as good.  A good way to see what has happened is to look at the sampling points:

Get the sampling points for numerical integration of f.

Block[{c = 0},rs = Reap[Sow[NIntegrate[f[x],{x,0,10}, EvaluationMonitor:> Sow[{x,++c}]]]][[2,1]]; c]

Get the sampling points for numerical integration of the black box version of f.

Block[{c = 0},rb = Reap[Sow[NIntegrate[blackbox[x],{x,0,10}, EvaluationMonitor:> Sow[{x,++c}]]]][[2,1]]; c]

Show the sample points for the two functions.

ListPlot[{rs, rb}]

In Mathematica 6, NIntegrate had a complete overhaul, including many new methods.  You can find a detailed description of methods and how to use them in the tutorial on  Advanced Numerical Integration in Mathematica.

Here is a final interesting example shown the sampling points used for a 3 dimensional integral.

points=Reap[NIntegrate[Boole[(x^2)+(((y-1))^2)<1]+Boole[(x^2)+(((z-1))^2)<1],{x,0,1},{y,0,2},{z,0,2},Method->{"GlobalAdaptive",Method->{"TrapezoidalRule","Points"->3},"SingularityHandler"->None},PrecisionGoal->4,EvaluationMonitor:>Sow[{x,y,z}]]][[2,1]];

Graphics3D[Point[points],Axes->True,AxesLabel->{"x","y","z"}]

Optimization

Unconstrained optimization

The functions FindMinimum and FindMaximum for doing unconstrained optimization provide a good example for how to control tolerances and computation precision with the AccuracyGoal, PrecisionGoal and WorkingPrecision options.

Start an unconstrained maximization from the point 0.

sc1 = FindMaximum[Sin[x] - Cos[x], {x,0}]

Starting at a different point leads to trouble.

sc2 = FindMaximum[Sin[x] - Cos[x], {x,-.56}]

The minimum is (Sqrt[2]) that occurs at x = ((3π)/4).

Check the actual error of the results.

{sc1[[1]] - (Sqrt[2]), (x /. sc1[[2]]) - 3 π/4}

{sc2[[1]] - (Sqrt[2]), (x /. sc2[[2]]) - 3 π/4}

So why has the error control run amok in the second case?

Looking at plots of the function at appropriate scales provides the answer:

s1 =Reap[FindMaximum[Sin[x] - Cos[x], {x,0}, EvaluationMonitor:>Sow[{x, Sin[x] - Cos[x] - Sqrt[2]}]]][[2,1]];
s2 =Reap[FindMaximum[Sin[x] - Cos[x], {x,-.56},EvaluationMonitor:>Sow[{x, Sin[x] - Cos[x] - Sqrt[2]}]]][[2,1]];

Manipulate[Show[Block[{ϵ = 10.^k, $DisplayFunction = Identity, mx = 3 π/4},{Plot[Sin[x] - Cos[x] - Sqrt[2], {x,mx -ϵ,mx+ϵ}, PlotLabel->"Width "<>ToString[InputForm[ϵ]],PlotPoints->1000, Method->{Compiled->False}],Plot[Block[{xp = SetPrecision[x, 1000]},Sin[xp] - Cos[xp] - Sqrt[2]], {x, mx -ϵ,mx+ϵ},PlotStyle->RGBColor[0,0,1], WorkingPrecision->1000, Method->{Compiled->False}]}],Graphics[{PointSize[0.02],{Green, Point[s1]},{Red, Point[s2]}}]],{{k,-6},-10,0,1}]

The fix for the problem is to allow higher computational precision:

FindMaximum[Sin[x] - Cos[x], {x, RandomReal[]},WorkingPrecision->20, AccuracyGoal->MachinePrecision/2, PrecisionGoal->MachinePrecision/2]

Constrained optimization

Mathematica now has a lot of methods for doing constrained optimization, from linear programming to nonlinear local constrained optimization to global optimization methods.  A good place to read in detail about this in the tutorial on  Constrained Optimization.

One issue when doing constrained optimization is to consider whether you want  a local or global minimum.  As an example, consider the following bumpy function.

Define a function s and show its surface over the unit circle in the plane.

Clear[a,b,c,d,e,f,g,h,s];
{a,b,c,d,e,f,g,h} = {5,0,7,1,4,3,9,2};
s[x_, y_] = Sin[a x^2 + b] + Cos[c y^2 + d] + Sin[e (x+ y)^2 + f] + Cos[g (x-  y)^2 + h];

p = Plot3D[s[x,y],{x,-1,1},{y,-1,1},RegionFunction->Function[{x,y},x^2 + y^2 ≤ 1], MeshFunctions->{#3&}]

FindMinimum finds a local constrained minimum if you include constraints.

Find a local constrained minimum of the function s.

{lm, lmxy} = FindMinimum[{s[x,y],x^2 + y^2 ≤ 1},{x, y}]

Show[p, Graphics3D[{Red, PointSize[0.05],Point[{x, y, lm} /. lmxy]}]]

This is a local minimum, but clearly not the global one.

NMinimize will work harder to attempt to find a global minimum.

Use the RandomSearch method of NMinimize to search for a global minimum.

{gm, gmxy} = NMinimize[{s[x,y],x^2 + y^2 ≤ 1},{x,y}, Method->"RandomSearch"]

Show[p, Graphics3D[{Red, PointSize[0.05],Point[{x, y, lm} /. lmxy]}],Graphics3D[{Green, PointSize[0.05],Point[{x, y, gm} /. gmxy]}]]

Rootfinding

Numerical rootfinding functions in Mathematica are NRoots and FindRoot.

I'll use rootfinding to demonstrate how different numerical functions can be used together effectively to get a result.

As an example, I'll show you how you can use FindRoot and NDSolve to implement a shooting method for boundary value problems.  The idea of a shooting method is to find an initial condition such that the solution passes through a specified boundary value.

Consider the BVP  y''(x) + sin(y(x)) == 0, y(0)==1, y(10) = 0

The first stage in setting up a shooting method is to define a function for FindRoot to work with.  In general, when you want to use one command inside of another, it is a really good idea to define a function that does it.  If nothing else, you can plot the function to see that it is working the way you expect.

Define a function of the derivative of y at x = 0 that gives the value of the solution y at x = 10. The function is defined to only evaluate when yp0 is a number since symbolic evaluation would be meaningless (and just plain wrong too)

Clear[x, y, t];
(Subscript[y, 10])[yp0_?NumberQ]:=First[y[10]/.NDSolve[{y''[x]+Sin[y[x]]==0,y[0]==1,y'[0]==yp0},y,{x,10}]]

Plot the function.

Plot[(Subscript[y, 10])[yp0], {yp0, -2,2}]

The function has multiple roots.  Each one corresponds to a different solution.

Find a root of the function.

Timing[r = FindRoot[(Subscript[y, 10])[yp0], {yp0,-1.75}]]

The default is to use Newton's method which has worked, but perhaps is not as fast as possible since the derivative needs to be evaluated using finite differences.  It is possible to define what is effectively a symbolic derivative by finding the ODE giving the sensitivity of the solution to the parameter y'(0).  An easier thing to do, however is to use a derivative free method.

Use Brent's method Find a root of the function.

Timing[r = FindRoot[(Subscript[y, 10])[yp0], {yp0, -2,-1.75},Method->"Brent"]]

sol = First[NDSolve[{y''[x]+Sin[y[x]]==0,y[0]==1,y'[0]==yp0/. r},y,{x,10}]];
Plot[Evaluate[y[x] /. sol],{x,0,10}, PlotRange->All]

In version 6, the shooting method is built into NDSolve:

bsol =First[NDSolve[{y''[x]+Sin[y[x]]==0,y[0]==1,y[10]==0},y,{x,10}]]

Plot[Evaluate[y[x] /. bsol],{x,0,10}, PlotRange->All]

Summary

I have chosen examples that show how numerical computation works in each of the (numerical) fields shown in the hierarchy above.  This is indeed a very quick overview of numerical computation in Mathematica, but should be enough to give you a fell of how the system was thought out as a whole and different tasks can be combined to get otherwise difficult results.  
Some of the important points I have illustrated are:

In operations based on input values or data, Mathematica automatically returns a result based on the precision of the inputs.

Functions are task or problem based rather than algorithm based.

Functions try to pick a good algorithm based on the problem given.

Method options give you detailed control over what algorithm is used and its parameters.

Output from one function can easily be used in other functions.