Introduction to Numerics

Introduction

Abstract

If all problems could be solved analytically, there would be no need for numerical approximations.  However, many important problems have no analytic solution, take too long to compute one, or are posed with inexact input.   The numerical functions in Mathematica offer a variety of sophisticated numerical algorithms that allow you to continue from symbolic calculations.

This notebook will give a brief overview of some of the numerical functions in Mathematica and show some of their typical uses.

This tutorial is written for Version 3.0 Mathematica.

Initialization

\(SetOptions[Plot, \ Frame -> True, Axes -> False, 
	GridLines -> None, 
      
	PlotStyle -> {Hue[0], RGBColor[0, 0, 1]}, ImageSize -> {400, 280}]; 
    SetOptions[ParametricPlot, \ Frame -> False, Axes -> True, 
	
      GridLines -> None, 
	PlotStyle -> {Hue[0], RGBColor[0, 0, 1]}, 
      ImageSize -> {400, 280}]; 

    SetOptions[ListPlot, Frame -> True, PlotRange -> All, PlotJoined -> True, 
      GridLines -> Automatic, ImageSize -> {400, 300}]; 

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

Numerical Expressions

Mathematica works with many different types of numerical expressions,  numbers and formulae that contain numbers.  This is one of its strengths.

Numbers

Two basic types of numbers, exact and approximate.

For example,

Sqrt[2]

is an exact quantity representing the positive root of (x^2) == 2.

When exact quantities are used in calculations, the results are exact:

(Sqrt[2])^2

Log[Sqrt[2]] - Log[2]/2

Numbers can typically be represented by decimals.  When the number of digits is restricted, we have an approximation.  For example,

approx = N[Sqrt[2]]

is an inexact approximation to the positive root of (x^2) == 2.    N calculates the numerical approximation.

When inexact numbers are used in calculations, the result are inexact:

approx^2

(Note the decimal point)

Log[approx] - Log[2]/2

These approximate numbers are machine numbers.  Mathematica can generate arbitrary precision approximations to numeric quantities.  The precision refers to the number of correct digits in the approximation.   Such numbers are often called bignums.

For example,

approx32 = N[Sqrt[2],32]

has 32 correct digits.

Now suppose we do a calculation, as before

approx32^2

It is obvious here that the result is correct, in other cases, you may want to know how many digits are likely to be correct.  Mathematica has a command, Precision, that gives you this information.

Precision[%]

Formulae

Mathematica  can also work with a formula,  a compound mathematical expression.

f = 
      1/((x - 0.2)^2 + 0.01) + 1/((x - 0.8)^2 + 0.04)

There are lots of things we can do with this formula,   plot it.

Plot[f, {x, -1, 2}, Frame -> True];

Differentiate it.

D[f, x]

Plot the derivative.

Plot[%, {x, -1, 2}, PlotRange -> All];

In the original picture we could see a minimum,   let's search for the minimum near 0.5.

FindMinimum[f,{x,0.5}]

A list results,  the first element is the minimum value,  the second the value of x.  We substitute in as follows.

f /. Last[%]

Something else we can do is quadrature...  area under the curve.

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

Summary

This is a basic part of Mathematica numerics,  so to summarize.

Finding Roots

FindRoot is the basic function for finding roots.

Computing the square root of 2, provides a good example.  Of course Mathematica has it's own square root function,
but this makes a good example.

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

Making a plot often helps see what we are doing.

Plot[(x^2), {x, 0, 3}, Prolog -> Line[{{0, 2}, {3, 2}}]];

The command above used Newton's method.  This is the best method to use when the derivative of the function can be computed directly.  Sometimes this is not possible, in which case it is best to use the secant method (Brent's method in 1D if the root is bracketed).  This is done by specifying two starting points for the algorithm.

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

The root is reasonably accurate:

x^2 /. %

Suppose you wanted to know how many digits of the root were correct.  This is easy to do in Mathematica by using bignums.  In FindRoot (and most other numerical functions as well), the precision of numbers used is controlled by the option WorkingPrecision.  For example, to work with 32 digits, you can use

FindRoot[x^2 == 2,{x,0,1},WorkingPrecision->32]

two = x^2 /. %

This does not have 32 digits correct:

two-2

But Mathematica still claims that the precision of "two" is 32.

Precision[two]

The discrepency is a part of the nature of numerical analysis.  In doing any numerical procedure, there are typically two types of errors.  Arithmetic (roundoff) error and truncation (algorithm) error.  The Mathematica command Precision keeps track only of the arithmetic error.  The truncation error is controlled by the algorithm.  In this case, Mathematica is using the secant method, which iteratively improves the result.  Mathematica must decide at some point when to stop the iterations.  This is controlled by the option AccuracyGoal.

The default AccuracyGoal->Automatic tells Mathematica to stop the iterations when the value of the function is zero accurate to WorkingPrecision-10 digits.  Why not go all the way to WorkingPrecision digits?  This may not always be possible arithmetically, as roundoff error may exceed the truncation error.  You can see that the result above  had the function value (10^(-26)), which is less than the stopping criterion (10^(-22)).

FindRoot will also do multidimensional root finding.  For example, consider the problem of finding the value of b for which the curve y = (b^x) is tangent to the line y = x.

Define the objective function.   (b^x) == xand D[ (b^x), x] == 1.     Use dot product to combine.

func  = {b^x-x,D[b^x,x]-1};

So find root of the dot product.

func . func

One way to obtain starting values is to use a plot to estimate where the root might be:

ContourPlot[func.func,{x,.1,4},{b,.1,4},
	PlotPoints->30,
	Contours->{0.1,1,2},
	PlotRange->{0,2}];

It looks like a good bet to try a value somewhere in the innermost contour.

FindRoot[func,{x,2},{b,1.5}]

A plot also verifies that we have the right solution:

Plot[{t,b^t /. %},{t,0,4}];

It would be possible to use FindMinimum to find this solution.

Linear Algebra

The symbolic/numeric hybrid approach applies also to linear algebra.  For large matrices,  symbolic calculations are generally impractical.  In this section we will look at some examples of numerical computation with matrices in Mathematica.

When we think about matrices we should bear in mind that many techniques can be generalized to more general tensors.  Tensors can be described with n-indices,  rather than the two indices that are used for matrices.   The Mathematica implementation is very general.

Constructing Matrices

Table is a basic constructor for matrices.     Here we get a 4×4 matrix of random values.

M = Table[Random[],{4},{4}]

Matrices are formed with Mathematica lists;  this means they work nicely with other parts of the system.  By default they do not format like matrices,  this can be done with MatrixForm.

M // MatrixForm

Looking at the prettified form can obscure some important details,  so generally it would not be good to wrap result from Table in MatrixForm,  before assigning to a variable.

We can just type the list notation in directly.

V = {1, 2, 3, 4}

Alternatively we can use Table to make a vector (a list of elements).

V = Table[i,{i,4}]

Here we make a rank-3 tensor ( a list of lists of lists).

T = Table[Random[],{4},{3},{3}]

There are many useful commands for getting information about matrices,  such as the dimensions and rank.

Dimensions[ T]

TensorRank[M]

Another way to generate a matrix, sometimes good for a symbolic matrix, is with Array.

Array[Subscript[m, SequenceForm[##]]&, {3, 3}] // MatrixForm

Array[Subscript[m, SequenceForm[##]]&, {3, 5, 3}] // MatrixForm

We can also generate diagonal and identity matrices.

IdentityMatrix[5] // MatrixForm

DiagonalMatrix takes a vector of the diagonal terms,  (here they are symbolic).

DiagonalMatrix[{a, b, c, d}] // MatrixForm

Accessing Elements

It is very important to be able to extract parts of a matrix.  This can be done with the standard Mathematica functions Part, Take, and Drop.

First we generate, and then print, a matrix for demonstration.

mat = Array[Subscript[m, SequenceForm[##]]&, {5, 5}];

mat // MatrixForm

Part

If we want to extract an element we use Part, entered here with a shortcut notation,

mat[[1, 2]]

alternatively we can type in the full name  Part,

Part[mat,  1, 2]

These two ways are completely identical.

This is a matrix so when we give Part two integer indices we extract an element.   If we gave one index we would get a row.

mat[[3]]

We can also use Part to pick out particular rows,  for example the (2^nd) and (3^rd) rows.

mat[[ { 2, 3}]]

This is really very powerful,  to get a particular column  we use a list to select all rows and an index for the column.  Here we get the (3^rd)column.

mat[[ Range[1, 5],  3]]

This works because we are taking rows 1 through 5,  ie all rows.

Extracting a Sub-matrix

To extract a sub-matrix  we can also use Part.

The sub-matrix from (Subscript[m, 21]) to (Subscript[m, 52]).

Part[ mat,  Range[2, 5], Range[1, 2]] // MatrixForm

We can use All to describe all the rows.  Negative indices count from the end.

This returns the matrix with the first and last columns dropped.

Part[ mat,  Range[1, 5], Range[ 2, 4]] // MatrixForm

We can also give a step size.

Part[mat, Range[1, 5, 2]] // MatrixForm

Take the first to the last column in steps of 2.

Part[mat, Range[1, 5], Range[1, 5, 2]] // MatrixForm

Tensors

All of this works quite properly with tensors,  not just with matrices.

mat = Array[Subscript[m, SequenceForm[##]]&, {2, 4, 2}];

mat // MatrixForm

This returns the vector found at position (Subscript[m, 32]).

Part[mat,  2, 2]

The sub-tensor from (Subscript[m, 221])to (Subscript[m, 223]).

Part[mat, {2}, {2, 3}, {1, 2}] // MatrixForm

Operations

Basic Arithmetic

The standard Mathematica operations such as Plus and Times are scalar operations rather than  the matrix operations.

mat = Array[Subscript[m, SequenceForm[##]]&, {2, 2}];

mat1 = Array[Subscript[n, SequenceForm[##]]&, {2, 2}];

mat // MatrixForm

If we multiply by a scalar this is applied to all elements.

10 mat // MatrixForm

If we multiply two matrices,  we get element by element multiplication.

mat mat1 // MatrixForm

To get matrix multiplication we use the Dot operation,  typically entered with a '.' .

mat  .  mat1 // MatrixForm

This of course works for general tensors.   For matrix multiplication of a rank-m and a rank-n tensor the result will be a rank-(m+n-2) tensor.    Consequently for matrix multiplication of two vectors the result is a scalar.

{a, b, c} . {x, y, z}

Matrix Functions

Starting with a random number matrix.

M = Table[Random[],{4},{4}]

If we want to invert this matrix,  since it was constructed randomly, we do not know a priori whether it is invertible or not.  One way to check is to see if it has  an empty null space.

NullSpace[M]

This means that we should be able to find its inverse.

inv = Inverse[M]

We can check on the inverse by taking the matrix product.

M.inv

The Chop command is useful to get rid of things that are near, but not quite zero:

Chop[%] // MatrixForm

Here we have a random vector b.

b = Table[Random[],{4}]

Suppose you want to solve the system M x ==  b for the vector x given the vector b.  One way is to use the inverse:

x1 = inv.b

Of course page 1 of any numerics book tells you not do this.

the correct way, both in terms of efficiency for large matrices and in terms of numerical stability,  is to solve the system by elimination.  This is done in Mathematica using LinearSolve

LinearSolve[M,b]

If you want to find the solution for many different vectors b, you can use the commands LUDecomposition and LUBackSubstitution as an alternative to finding the inverse:

dec = LUDecomposition[M];
LUBackSubstitution[dec,b]

The eigenvectors and eigenvalues are easily found:

Eigenvalues[M]

Eigenvectors[M]

Numeric Functions

A basic property of many Mathematica functions is that they are Listable,  this means they will thread over the elements of a matrix or tensor.

mat = Array[Subscript[m, SequenceForm[##]]&, {4, 4}];

The Sin function is Listable,  this means that it will thread down over every list (row and column) in a matrix.

Sin[mat] // MatrixForm

You can make your own functions Listable.

SetAttributes[ fun,  Listable];

fun[mat] // MatrixForm

Here you still need to say how fun should behave for each element.

It is very advantageous to use listability to operate on the elements of a matrix.   For numerical matrices these computations can be done in a very efficient manner,   one can consider this as a form of parallel or vectorised computation.

Setting Elements

We can use Part on the left-hand side of a set to change elements.

mat = Array[Subscript[m, SequenceForm[##]]&, {5, 5}];

mat // MatrixForm

Set the element 3,3 to be x.

mat[[3, 3]] = x;

mat // MatrixForm

Here we set a row to have a particular value...

mat[[2]] = {1, 2, 3, 4, 5};

mat // MatrixForm

Here we set a column to have a particular value...

mat[[Range[1, 5], 2]] = {-1, -2, -3, -4, -5};

mat // MatrixForm

Vectors

A final comment about vectors in Mathematica.

Matrices in Mathematica are represented by a list of lists of elements.  Vectors are represented by a list of elements.

Sometimes a basic question comes up.

'How do I represent a row and a column vector in Mathematica?'

Answer:   'there are no such things'.

You can have a matrix with one row or one column...

row = {{x, y, z}}; \(
\)col = {{x}, {y}, {z}};

row // MatrixForm

col // MatrixForm

but these are matrices...  if I multiply the row by the column I do not get a scalar,  I get a matrix.

row . col

Really if I take the vector product of two vectors I get a scalar.  Of course you need to distinguish between them.   As Mathematica does,  it is not constructive to do otherwise.

If you are unconvinced ponder this.   The points in an n-dimensional space can be represented with a length-n vector.   Is this a row or a column vector?   If I have a plane in this space,  is the normal to this a row or a column vector?   The questions are meaningless.   A length-n vector is just a collection (an n-tuple) of elements.

Differential equations

The command NDSolve will solve (systems of) differential equations.  For example, Duffings equations has no analytic solutions, but can easily be solved with NDSolve:

sol = NDSolve[{
	x''[t] + .15 x'[t] + x[t] - x[t]^3 == .3 Cos[t],
	x'[0] == x[0] == 0},
	x,{t,0,10}]

The InterpolatingFunction object returned can be treated like any Mathematica function that evaluates to a number.   Here we make a phase space plot, easy in Mathematica since we just differentiate the result of NDSolve.   Good example of symbolic/numeric powers.

ParametricPlot[Evaluate[{x[t],x'[t]} /. sol],{t,0,10}];

The methods used in NDSolve include a solver for "stiff" systems of equations.  A system of equations is called stiff when it has components that vary on extremely different time scales.  For example,

sol = NDSolve[{
	x'[t] == -1000 x[t] + Cos[t],
	y'[t] == -y[t] + t^2, 
	x[0]==1, y[0] == 0},
	{x,y},{t,0,2}];

Plot[Evaluate[{x[t],y[t]} /. sol[[1]]],{t,0,2},
	PlotRange->All];

The first component of the system varies much more rapidly.

Plot[Evaluate[x[t] /. sol[[1]]],{t,0,.005}];

Data Handling

Mathematica handles vectors as lists and matrices as lists of lists.  We can work with data using lists.  There are many commands to work with data.  Brief description here.

Importing Data

ReadList is probably the main function for importing data.

Data Analysis

Linear fitting can be done with Fit.

data = 
      Table[{x, Exp[-x] Sin[x] + (Random[] - 0.5)/20. }, {x, 0.1, 
          10, 0.05}];

p1 = ListPlot[data];

fun = Fit[data, {1, x, x^2, x^3}, x]

p2 = Plot[fun, {x, 0, 10}, PlotRange -> All];

Show[p2, p1];

Linear fit is not very satisfactory.  Lets try a non-linear fit,  from the appropriate package.

<< Statistics`NonlinearFit`

Plot[ Exp[-x]  -  Exp[-2 x], {x, 0, 10}];

b

NonlinearFit[data,   Exp[-b1 x] -  c1 Exp[-d1 x], 
      x, {{b1, 1}, {c1, 1}, {d1, 2}}]

p2 = Plot[%, {x, 0, 10}];

Show[p1, p2];

Statistics

There are many standard packages that carry out a wide variety of statistical computations.  In addition there are application packages to provide functions in TimeSeries or Wavelets.   Really this is outside of the scope of this tutorial.   Look in the on-line documention.

Quadrature and Summation

Quadrature

The command NIntegrate can approximate numerical integrals accurately even over infinite intervals.

NIntegrate[Exp[-x^2/2]/Sqrt[2 Pi],{x,-Infinity,Infinity}]

NIntegrate assumes its integrand is smooth.

Clear[f];
f[x_] := Which[
	x < 0,	0,
	x < 1, 	x^2,
	x < 2,	2x,
	x < 3,	4,
	x >= 3,	0]

NIntegrate[f[x],{x,-3,6}]

This does quite poorly bacause the assumption of smoothness is violated.  You can tell NIntegrate to treat different regions separately by including the points where there might be discontinuities.

NIntegrate[f[x],{x,-3,0,1,2,3,6}]

NIntegrate can handle singularities at the endpoints of its region (or in the middle if the points are specified as above).  For example,

NIntegrate[1/Sqrt[Abs[x]],{x,-1,0,1}]

Sums and Products

There are many ways to compute sums in Mathematica.  I will use a simple example to demonstrate the difference between these.

Timing[N[Sum[Exp[-x]Sin[x],{x,1,10000}]]]

N[Sum[]] works first by summing the sum (symbolically or numerically depending on the arguments) and then applying N to the result.

Timing[Sum[N[Exp[-x]Sin[x]],{x,1,10000}]]

Sum[N[],] adds up the numerical values of the summands.

Timing[NSum[Exp[-x]Sin[x],{x,1,10000}]]

NSum works by adding up the first few terms, then using mathematical approximations to estimate the size of the tail.  Even though it was much faster, its result is just as good in this case:

%[[2]]-%%[[2]]

Because NSum approximates the tail, it can do infinite sums as well:

NSum[1/(n(n+1)(n+2)),{n,1,Infinity}]

In cases where the Sum cannot be evaluated, this is equivalent to N[Sum[]]

N[Sum[1/(n(n+1)(n+2)),{n,1,Infinity}]]

However, if you know what you want is a numerical approximation to the sum, NSum has the advantage that it takes options that control the way the sum is approximated.

NSum uses two different methods to approximate the tail.  One is Euler-MacLaurin summation, where the tail is approximated by integration.  This was the method used in the above sum.  The other method uses Wynn's epsilon algorithm (as in the Mathematica command SequenceLimit) to approximate the limit of the sequence of partial sums.  This method is used in sums like:

NSum[1/2^n,{n,1,Infinity}]-1

Mathematica uses heuristics to decide which method is best to use.  You can control the method using the option Method

NSum[1/2^n,{n,1,Infinity},Method->Integrate]-1

In this case, Integrate was not a great choice.

NProduct works similarily to NSum, but it computes the exponential of the sum of the logarithms of the terms instead.

NProduct[1+2^-n,{n,1,Infinity}]

Summary

These are some of the key elements of Mathematica numerics.