Introduction to Numerical Computation in Mathematica

Rob Knapp

Wolfram Research, Inc.

Abstract

Mathematica is, in part, a very general system for numerical computation. Functionality for numerical computation range from simple arithmetic to sophisticated optimization and integration functions.  A strength of Mathematica is that a common set of principles underlies the design of the functions. These common principles will be presented and then demonstrated in a quick overview of the numerical functionality.

Initialization

Off[General::spell];
Off[General::spell1];

SeedRandom[1234];
data  =  
        Block[{a  =  2,  b  =  1, x},  
          Table[{x  =  4 Random[],  
              Random[]  +  
                NIntegrate[
                  Exp[a Sin[b s]] ,  {s, 0, x}]},  {50}]];

Introduction

To start off, I'll give a few quick examples of numerical computations done with Mathematica.  Then I will introduce some common design principles underlying all of the numerical functionality in Mathematica and give more extensive examples exploring these principles.

(Sqrt[2])

For the first few examples, I'll show you several different ways of computing (Sqrt[2]) numerically.

N gives the numerical value with machine number arithmetic.

N[(Sqrt[2])]

N[e, p] gives the numerical value using arbitrary precision software arithmetic.

N[(Sqrt[2]),  100]

Find (Sqrt[2]) as a root of a simple equation.  FindRoot uses an iterative search from a given starting point.

FindRoot[x^2  ==  2,  {x, 1}]

If you just want the number, you can use ReplaceAll (/.)

x  /.  %

Find (Sqrt[2]) as a root of a simple equation to higher precision.

FindRoot[x^2  ==  2,  {x, 1. },  
      WorkingPrecision -> 100]

Find (Sqrt[2]) as an largest eigenvalue of a matrix.

Eigenvalues[{{1. , 1. }, {1. , -1. }}, 1]

Find (Sqrt[2]) as an largest eigenvalue of a matrix of higher precision numbers.

Eigenvalues[{{1.`100, 1.`100}, {1.`100, -1.`100}}, 1]

Find (Sqrt[2])the first time where the solution of a differential equation is zero.

sol  =  
      NDSolve[{x''[t]  +  (1/2) (π^2) x[t]  ==  0,  
          x[0]  ==  0,  x'[0]  ==  1},  
        x,  {t, 0, ∞},  
        Method -> {"EventLocator",  "Event" -> 
              x[t]}]

Plot[Evaluate[First[x[t]  /.  sol]], {t, 0, (Sqrt[2])}];

Find (Sqrt[2]) from twice the value of t where the solution just found has its maximum value.

f[t_?NumericQ]  :=  
        First[x[t]  /.  sol];
FindMaximum[f[t],  {t, 0}]

2*(t  /.  %[[2]])

Numerical function design principles

Algorithm selection

I could spend the entire time on this, and probably for just one type of problem.  But to keep moving along, I will give just a couple of relatively simple examples.

Find a solution to an unconstrained minimization problem.

FindMinimum[(x -  1)^2  +  
        100 (y  -  Sin[1.5 Pi x])^2,  {{x, -1}, {y, 
          1}}]

You might be wondering how FindMinimum arrived at this solution.  For the commands that work with function evaluations, there are options EvaluationMonitor and StepMonitor that allow you to keep track of where the algorithm did sampling.

See where the algorithm was at each step.  Reap and Sow are functions that can be used to collect data in a convenient way.

Reap[
      FindMinimum[(x -  1)^2  +  
          100 (y  -  Sin[1.5 Pi x])^2,  {{x, -1}, {y, 1}},  
        StepMonitor :>  Sow[{x, y}]]]

This shows the "step" points.

lp  =  ListPlot[%[[2, 1]]];

Show the step points along with a contour plot of the function.

Show[
        Block[{$DisplayFunction  =  Identity},  { 
            ContourPlot[(x -  1)^2  +  
                100 (y  -  Sin[1.5 Pi x])^2,  {x, -1, 
                1}, {y, -1, 1},  ContourShading -> True,  
              ColorFunction -> Hue], lp}]];

There is a package that does the same thing, but making a fancier plot:

Needs["Optimization`UnconstrainedProblems`"]

FindMinimumPlot[(x -  1)^2  +  
        100 (y  -  Sin[1.5 Pi x])^2,  {{x, -1}, {y, 
          1}}]

Solve the problem using Newton's method.

FindMinimumPlot[(x -  1)^2  +  
        100 (y  -  Sin[1.5 Pi x])^2,  {{x, -1}, {y, 1}},  
      Method -> "Newton"]

Solve the problem using a Quasi-Newton method.

FindMinimumPlot[(x -  1)^2  +  
        100 (y  -  Sin[1.5 Pi x])^2,  {{x, -1}, {y, 1}},  
      Method -> "QuasiNewton"]

Solve the problem using a derivative-free method.

FindMinimumPlot[(x -  1)^2  +  
        100 (y  -  Sin[1.5 Pi x])^2,  {{x, -1, - .9}, {y, 
          1,  .9}},  Method -> "PrincipalAxis"]

Numerically integrate f(x) = (e^(cos(4 x)))

{v,  {e}}  =  
      Reap[NIntegrate[Exp[Cos[4 x]],  {x, -1, 1},  
          EvaluationMonitor :> Sow[x]]]

Show where the function was evaluated.

ListPlot[e];

Numerically integrate f(x) = (e^(cos(4 x)))with the symbolic form hidden.

f[x_?NumberQ]  :=  Exp[Cos[4 x]];
{v,  {e}}  =  
        Reap[NIntegrate[f[x],  {x, -1, 1},  
            EvaluationMonitor :> 
              Sow[x]]];
v

Show where the function was evaluated.

ListPlot[e];

Propagation of precision and accuracy

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

Precision and Accuracy

The functions Accuracy and Precision allow you monitor this.  They give you information about the absolute and relative size the interval within which the value is not known with certainty.  You can roughly think of precision as the number of significant digits and accuracy as the number of significant digits to the right of the decimal point.

Machine numbers have precision MachinePrecision

Precision[1.0]

Numbers that are specified exactly have infinite precision.

Precision[(Sqrt[2])]

When given exact input, N works adaptively to obtain the requested precision.

Precision[N[(Sqrt[2]),  100]]

In most instances, Precision gives you a sufficient measure of how well you know a particular number. However, for numbers that could be zero, the whole concept of relative error breaks down.

(Sqrt[2]) + (Sqrt[3]) - (Sqrt[5 + 2 (Sqrt[6])]) is an obscure way to express 0.

nz  =  N[(Sqrt[2]) + (Sqrt[3]) - (Sqrt[5 + 2 (Sqrt[6])]),  20]

No significant digits are known (because it doesn't have any!)

Precision[nz]

However, you can still get a measure of the error using Accuracy.

Accuracy[nz]

Propagation of error

In most computations done with approximate numbers, error propagates with each operation.  In many cases, errors cancel themselves out and do not grow excessively.  However, there are cases where the error does become problematic.  I'll show examples from these cases, not because they are the most common, but because they highlight these aspects of Mathematica.

With machine number computations,  all numbers have the same precision, MachinePrecision, and to keep track of potential precision lost, it may be necessary to use error bounds or other techniques from numerical analysis for a particular algorithm of computation.

Do an iteration of the discrete logistic map with machine number arithmetic.

ml  =  NestList[4 # (1  -  #) &,   .1,  100]

Each number has the same precision.

Map[Precision,  ml]

Mathematica's arbitrary precision arithmetic is implemented in software, so it is possible to keep track of how much error might be made with each arithmetic operation.  Mathematica does this using a method called significance arithmetic.

Do an iteration of the discrete logistic map with arbitrary precision arithmetic.

al  =  
      NestList[4 # (1  -  #) &,  N[1/10,  50],  100]

The precision decreases with each iteration until there is none

Map[Precision, al]

Plot the precision shown by the two computations.

Show[
        Block[{$DisplayFunction  =  Identity},  {ListPlot[
              Map[Precision,  ml],  PlotStyle -> Red],  
            ListPlot[Map[Precision,  al]]}]];

Show the values in the two lists where they actually cross.

al[[60]]

ml[[60]]

Significance arithmetic uses approximations to determine the propagation of error.  For numbers with even moderate precision, the approximations are quite accurate.  However, when the precision becomes very small, these approximations can break down.  Mathematica also has an interval arithmetic that is slower, but more rigorously correct, particularly as intervals get bigger.

Fixed precision

For some computations, it is appropriate to keep the precision fixed, as it is in machine hardware arithmetic.  

For example, Mathematica does most of its numerical linear algebra with fixed precision arithmetic because the most robust and efficient algorithms have been developed based on error bounds that depend on a fixed precision model.

Here is a Hilbert matrix.

TableForm[
      m  =  Table[
          1/(i  +  j  -  1),  {i,  10},  {j,  10}]]

Choose a right hand side for a linear system so the solution should have all components equal to 1.

b  =  m . Table[1, {10}]

With exact input, Mathematica gives an exact result.

LinearSolve[m, b]

If any of the inputs are at machine precision, the result is computed with machine arithmetic.

LinearSolve[m,  N[b]]

You can do computations of your own in fixed precision by setting the system parameters $MinPrecision and $MaxPrecision equal to each other.  To avoid potentially messing up other calculations, I recommend that you always do this within Block.

Do the logistic iteration with fixed 50 digit precision.

fl  =  
      Block[{$MinPrecision  =  50,  $MaxPrecision  =  
            50}, 
NestList[4 # (1  -  #) &,  
          1/10`50,  100]]

Show a plot of the differences.

Off[Graphics::gptn];
Show[
        Block[{$MinPrecision  =  50,  $MaxPrecision  =  
              50,  $DisplayFunction  =  
              Identity}, 
{ListPlot[
              Log[10,  Abs[fl  -  al]]],  
            ListPlot[Log[10,  Abs[fl  -  ml]],  
              PlotStyle -> Red]}]];

Intervals

innerval  =  N[Interval[1/10], 50]

fint  =  
      Block[{$MinPrecision  =  50,  $MaxPrecision  =  
            50}, 
NestList[4 # (1  -  #) &,  
          innerval,  100]]

Log[10. ,  Apply[#2  -  #1 &,  fint[[All, 1]], 1]]

Quality Control

For numerical algorithms that work with functions, it is not as straightforward to get a result with a certain precision or accuracy as it is in arithmetic.  These functions allow to specify desired tolerances in terms of goals, in particular PrecisionGoal and AccuracyGoal, representing relative and absolute tolerances.

Alas, it is not quite as simple as specifying the goal and being sure that you'll get a result correct to that tolerance.  First, numerical algorithms of this type typically work by sampling the function, so can miss features entirely., or, if the function itself is not computed accurately enough the algorithm may not work properly.   Second, error measure is based on estimation from information built up during the running of the algorithm and those estimates are not the same as the actual error.  Finally, different algorithms have different types of error control.  I'll show three types of error control, convergence estimates, global error control, and local error control.

Iterative convergence estimates

Many of the optimization functions work by iterations that eventually converge to the desired result.  Error estimates are typically based on the assumption (usually verified in some way) that the algorithm is converging and that the function is reasonably smooth.

Yet another way to compute (Sqrt[2])

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

Starting at a different point leads to trouble.

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

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?

A sequence of plots that explains the problem.

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

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

Global error estimates

When algorithms are able to keep track of global error estimates, it makes sense to do so.  In Mathematica, NIntegrate keeps track of the total (estimated) error for the entire integration problem and continues the algorithm until this is within tolerance.

Numerically integrate(1/(2 (Sqrt[ x]))) from 0 to 2 (yet another way to get (Sqrt[2]))

NIntegrate[(1/(2 (Sqrt[ x]))) , {x, 0, 2}]

%  -  Sqrt[2]

Numerically integrate(1/(2 (Sqrt[ x]))) from 0 to 2, trying to get a better value for (Sqrt[2])

NIntegrate[(1/(2 (Sqrt[ x]))) , {x, 0, 2},  
      PrecisionGoal -> 16]

%  -  (Sqrt[2])

If you think about it, its not surprising why this one ran into trouble: the desired precision is about the same as precision of the machine numbers being used to do the computation, leaving absolutely no room for error in the function values or arithmetic.

Numerically integrate(1/(2 (Sqrt[ x]))) from 0 to 2 using higher precision numbers, trying to get a better value for (Sqrt[2])

NIntegrate[(1/(2 (Sqrt[ x]))) , {x, 0, 2},  PrecisionGoal -> 16,  
      WorkingPrecision -> 24]

%  -  (Sqrt[2])

Even with relatively robust global error control, you can miss features that have small support.

TableForm[
      Table[{k  =  (10^i),  
          NIntegrate[
            (1/(2 (Sqrt[ x])))  +  (Sqrt[(k/π)]) 
                Exp[-k (((x -  (1/π)))^2)], {x, 0, 2}],  
          N[Integrate[
              (1/(2 (Sqrt[ x])))  +  (Sqrt[(k/π)]) 
                  Exp[-k (((x -  (1/π)))^2)], {x, 0, 2}]]}, {i, 1, 
          6}]]

Local error control

In some algorithms, it is more practical to control the error locally and hope that the local errors do not combine in a bad way to make the overall error too big.  NDSolve is an example of a function with local error control.  For an initial value problem, it starts at one end of the interval and takes a sequence of steps to get to the other end of the interval with each step constrained so that local error estimates are within tolerances.

Find the solution for the harmonic oscillator for a long time.

Timing[
      sol  =  NDSolve[{x''[t]  +  x[t]  ==  0,  
            x[0]  ==  1,  x'[0]  ==  0},  
          x,  {t, 0,  10000 π},  
          MaxSteps -> ∞]]

This shows the error at the end of each cycle.

ListPlot[
        Table[{t,  First[x[t]  /.  sol]  -  1},  {t,  0. ,  
            10000.  π,  2 π}]];

Increasing AccuracyGoal and PrecisionGoal takes longer ...

Timing[
      sol10  =  
        NDSolve[{x''[t]  +  x[t]  ==  0,  
            x[0]  ==  1,  x'[0]  ==  0},  
          x,  {t, 0,  10000 π},  MaxSteps -> ∞,  
          PrecisionGoal -> 10,  AccuracyGoal -> 10]]

... and it improves the solution.

ListPlot[
        Table[{t,  First[x[t]  /.  sol10]  -  1},  {t,  0. ,  
            10000.  π,  2 π}]];

Sometimes different methods will do better on problems with structure suited to that method.

Timing[
      solirk  =  
        NDSolve[{x''[t]  +  x[t]  ==  0,  
            x[0]  ==  1,  x'[0]  ==  0},  
          x,  {t, 0,  10000 π},  Method -> ImplicitRungeKutta,  
          MaxSteps -> ∞]]

This shows the error at the end of each cycle.

ListPlot[
        Table[{t,  First[x[t]  /.  solirk]  -  1},  {t,  0. ,  
            10000.  π,  2 π}]];

Access to data structures and lower level architecture

The numerical algorithms in Mathematica use a lot of different data structures.  It has many times been that case that a lower level data structure has proven very useful when made accessible through an integrated interface.  Two examples of this are PackedArrays and SparseArrays.

Packed Arrays

Packed arrays, introduced in Mathematica 4  are a data structure that we have been able to use to make many numerical functions much more efficient.  They are just an alternative representation of list that allow for faster access for operations and less storage space.

This is a list

r  =  Range[0,  20]

Actually, it is a packed array:

Developer`PackedArrayQ[r]

ByteCount[r]

A packed array can only have one machine native type (Integer, Real, or Complex) for all of its elements.

rr  =  1/N[r]

Developer`PackedArrayQ[rr]

r and N[r] use much less memory than rr

Map[ByteCount,  {r,  N[r],  rr}]

Mathematica generates packed arrays for larger lists when it is typically faster to have packed arrays rather then lists.

Generate a table of values of the sine function.  The Table command is automatically compiled to produce a packed array.

Timing[
      Developer`PackedArrayQ[
        Table[Sin[x],  {x,  0. ,  2.  Pi,  
            2.  Pi/1000000}]]]

Packed arrays propagate through arithmetic when possible.

Generate a table of values of the sine function by first generating the x values, then using vector arithmetic to get the sine values.

Timing[
      Developer`PackedArrayQ[x  =  N[Range[0, 1000000]]/(2.  Pi);  
        Sin[x]]]

Sparse Arrays

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]

Now lets make a big one:

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

Timing[f  =  LinearSolve[S]]

f[S . Range[1000]]

Combining commands

Something to keep in mind when combining commands is that in the inner functions, you may need to compute to somewhat stricter tolerances than your goal for the overall problem is.

Boundary value problems (FindRoot and NDSolve)

One way to find the solution to a boundary value problem is by shooting, where you try to find initial conditions will land you right on the target of the far boundary.  This example shows how you might do this.

Define a function of the derivative of y at x = 0 that gives the value of the solution y at x = 10.

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}]]

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.  
Note that I have used PatternTest (_?NumberQ) because the function would give nothing useful for symbolic yp0.

Plot the function.

Plot[(Subscript[y, 10])[(Subscript[y, 0])'],  {(Subscript[y, 0])',  -2, -1.75}];

Find a root of the function.

r  =  FindRoot[(Subscript[y, 10])[(Subscript[y, 0])'],  {(Subscript[y, 0])',  -1.75}]

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

Fitting data to an integrated function

In this example, I find a fit to data based on a model that is the integral of a function.

Show the data.

ListPlot[data];

The model function.

Clear[a, b];

model[a_?NumberQ,  b_?NumberQ][x_?NumberQ]  :=  
      NIntegrate[Exp[a Sin[b s]] ,  {s, 0, x}]

Finding the best fit.

Timing[
      fit  =  FindFit[data,  model[a, b][x], {a, b},  x,  
          PrecisionGoal -> 4,  AccuracyGoal -> 4]]

Show the model and the data together.

Show[
        Block[{$DisplayFunction  =  Identity}, {
Plot[
              Evaluate[model[a, b][x]  /.  fit], {x, 0, 
                4}], 
ListPlot[data]}]];

This worked, but it could have been done more efficiently.

model1[a_?NumberQ,  b_?NumberQ][x_?NumberQ]  :=  
      First[y[x]  /.  
          NDSolve[{y'[s]  ==  Exp[a Sin[b s]], 
              y[0]  ==  0},  y,  {s, 0, x}]]

Timing[
      fit  =  FindFit[data,  model1[a, b][x], {a, b},  x,  
          PrecisionGoal -> 4,  AccuracyGoal -> 4]]

Well, that's slower, not more efficient!  

The thing to observe is that with NDSolve, if you integrate once out to the maximum value of x needed, you can just use the InterpolatingFunction for subsequent evaluations of the model at different x for the same parameter values.  This can be done by rearranging rules.

maxx  =  Max[data[[All, 1]]];
model2[a_?NumberQ,  b_?NumberQ]  :=  
      First[y  /.  
          NDSolve[{y'[s]  ==  Exp[a Sin[b s]], 
              y[0]  ==  0},  y,  {s, 0, maxx}]]

model2[2.06,  1.02]

Timing[
      fit  =  FindFit[data,  model2[a, b][x], {a, b},  x,  
          PrecisionGoal -> 4,  AccuracyGoal -> 4]]

Even worse yet!

The whole point of using NDSolve was to take advantage of generating the InterpolatingFunction once and reusing it for each value of x for a given a,b value.  However, in the computation we just did, it is being recomputed for every a,b, and x.  It is easy to have Mathematica save the values computed for a given a,b:

model2[a_?NumberQ,  
        b_?NumberQ]  :=  (model2[a, b]  =  
          First[y  /.  
              NDSolve[{y'[s]  ==  Exp[a Sin[b s]], 
                  y[0]  ==  0},  y,  {s, 0, maxx}]])

Timing[
      fit  =  FindFit[data,  model2[a, b][x], {a, b},  x,  
          PrecisionGoal -> 4,  AccuracyGoal -> 4]]

In the examples I have done, I've tried to give you an introduction to numerical computing in Mathematica centered around some important design principles that underlay the numerical functions in Mathematica.