Introduction to Numerical Computing with Mathematica

Rob Knapp

Wolfram Research, Inc.

Initialization

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

Needs["GUIKit`"]

Needs["NumericalMath`EquationTrekker`"]

Needs["Optimization`UnconstrainedProblems`"]

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

The topic of numerical computing in Mathematica covers a lot of territory.  This presentation, being an introduction, does not purport to cover anywhere near everything involved, let alone a significant fraction.  Topics and examples have been chosen, however, to provide a broad introductory overview.

Even though I will be covering somewhat different topics, there are some common threads which run through much of what I will be demonstrating today.

Numbers and Arithmetic

These make the foundation upon which all numerical computations in Mathematica are built.

Types of numbers

Many systems provide only two or three types of numbers, typically machine sized integers, reals and complex numbers.   While this is simpler, such limited number systems do not have the flexibility of Mathematica, which supports integers, rationals, reals, and complex numbers of arbitrary size as well as exact representations of irrational numbers.

Exact numbers are integers, rationals, and representations for irrational numbers with no representation error.

Here is a large integer. Mathematica can represent integers with millions of digits.

1234567891011121314151617181920212223242526272829

Here is a rational number.

123456789101112131415/1617181920212223242526272829

Here is a representation for an irrational number.

irr  =  First[x  /.  Solve[x^2  ==  2, x]]

Mathematica has a function that allows you to determine whether an expression is a number or not.

A rational is stored as a number.

NumberQ[123456789101112131415/1617181920212223242526272829]

The irrational number, however, is a representation for a number as opposed to actually being a number.

NumberQ[irr]

It is important to be able to tell when an expression represents a number. This information can be obtained from Mathematica using the NumericQ  function.

This shows that irr, as assigned above, represents a number.

NumericQ[irr]

Exact representations of numbers are ideal in that you can combine them however you like and you get no error. However, since they are only representations, basic information, such as magnitude, sign, and so on, are not always apparent. With approximate numbers, such information is always obvious but with the cost of approximation error.

Here is an approximate number.

x = 1.23456789

As a default, Mathematica uses machine hardware numbers to store approximate numbers.  On the computer platforms Mathematica supports, these are IEEE 53 bit mantissa numbers.  The advantage of machine numbers is that they are efficient for both memory and speed.  However, because of the fixed size, there is a limited number of significant digits and size.

Here is an approximate number with too many digits specified to fit in a machine sized number.

y = 
      123456789.1011121314151617181920212223242526272829  +  
        8319837987398179379.1830813893809890891893 I

Here is an approximate number too large to fit in a machine sized number.

z = Exp[1000. ]

A simple speed comparison is useful to see the difference in efficiencies.

This gives the time it takes to square the machine number x a million times.

Timing[Do[x*x,  {10^6}];]

Numerical coercion

Mathematica allows coercion from exact numbers or numeric representations to approximate numbers through the command N.

This gives an approximate value for the square root of -2.

N[Sqrt[-2]]

This gives an approximate value for the factorial of 1000.

N[1000!]

The default is to use machine numbers when possible.  However, if you want to get a result with a particular precision, you can specify the precision in the second argument of N.

This gives an approximate value for the square root of -2 with 100 digits of precision.

N[Sqrt[-2],  100]

Significance arithmetic

A big advantage of Mathematica's software approximate numbers over machine numbers is that when arithmetic operations are done, not only is the result computed, but the precision and/or accuracy of the result is computed using an approximation known as significance arithmetic.

It is beyond the scope of this talk to describe the mathematics of significance arithmetic, but I will show a few examples to demonstrate how it works.  For more detailed information, I recommend reading the first part of Chapter 3 in The Mathematica Book.

The commands Precision and Accuracy give the significance of a number.  They have formal definitions in terms of error intervals they relate to.  It is often useful in terms of intuition, though not always strictly correct, to make an analogy between the precision and the number of significant (decimal) digits.

This gives an approximate number which is equal to one with 30 significant digits.

x  =  N[1,  30]

This takes a step of Newton's method for f(x) = (x^2)- 1 and checks the significance of the difference between the result and (Sqrt[2])

TableForm[{{x  =  x/2  +  1/x,  Precision[x]}, {(Sqrt[2])  -  x,  
          Precision[(Sqrt[2])  -  x]}}]

This takes a step of Newton's method for f(x) = (x^2)- 1 and checks the significance of the difference between the result and (Sqrt[2])

TableForm[{{x  =  x/2  +  1/x,  Precision[x]}, {(Sqrt[2])  -  x,  
          Precision[(Sqrt[2])  -  x]}}]

This takes a step of Newton's method for f(x) = (x^2)- 1 and checks the significance of the difference between the result and (Sqrt[2])

TableForm[{{x  =  x/2  +  1/x,  Precision[x]}, {(Sqrt[2])  -  x,  
          Precision[(Sqrt[2])  -  x]}}]

This takes a step of Newton's method for f(x) = (x^2)- 1 and checks the significance of the difference between the result and (Sqrt[2])

TableForm[{{x  =  x/2  +  1/x,  Precision[x]}, {(Sqrt[2])  -  x,  
          Precision[(Sqrt[2])  -  x]}}]

This takes a step of Newton's method for f(x) = (x^2)- 1 and checks the significance of the difference between the result and (Sqrt[2])

TableForm[{{x  =  x/2  +  1/x,  Precision[x]}, {(Sqrt[2])  -  x,  
          Precision[(Sqrt[2])  -  x]}}]

This takes a step of Newton's method for f(x) = (x^2)- 1 and checks the significance of the difference between the result and (Sqrt[2])

TableForm[{{x  =  x/2  +  1/x,  Precision[x]}, {(Sqrt[2])  -  x,  
          Precision[(Sqrt[2])  -  x]}}]

This takes a step of Newton's method for f(x) = (x^2)- 1 and checks the significance of the difference between the result and (Sqrt[2])

TableForm[{{x  =  x/2  +  1/x,  Precision[x]}, {(Sqrt[2])  -  x,  
          Precision[(Sqrt[2])  -  x]}}]

Significance arithmetic is very useful for showing when your results may have limited value.

This computes the sine of a large approximate number.

Sin[x^123]^2

This iterates the logistic map.

TableForm[
      NestList[(7/2) # (1  -  #) &,  1/x,  60]]

Vectors, Matrices, and Numerical Linear Algebra

Vectors and matrices are an important part of larger numerical computations.  Mathematica makes it easy to work with vectors and matrices and provides robust and fast numerical linear algebra functionality providing a second tier building block for numerical computations in Mathematica.

Vectors

A List is a representation for a vector.  In Mathematica, Lists are like other data structures but have special rules for arithmetic and some other functions.

Here is a list of numbers.

v  =  {1. ,  2,  3.   +  4.  I}

It is a vector.

VectorQ[v]

Vectors can be generated many ways but a very common and easy way of producing a vector is the Table command.

This makes a vector of squares of length 100.

tab  =  Table[N[i^2], {i,  250}]

You can get single or multiple parts of vectors with  Part and Take.

This gets the third element of the table.

tab[[3]]

This gets the first and last elements of the table.

tab[[{1, -1}]]

This gets every other element from the table.

Take[tab,  {1, -1, 2}]

Probably the most important property about lists is Listability.  Functions which have the Listable attribute act on list arguments elementwise.

This exponentiates with base 2 each element in the vector v.

2^v

This does a combination of listable operations.

Exp[2 v  -  {1, 2, 3}/v]  +  v v

Not only is it usually clearer and more compact notationally to use vector operations, but in general, it is much faster to do operations on vectors when possible.

This computes the square root of each element of the table using an explicit loop.

Table[Sqrt[tab[[i]]],  {i,  Length[tab]}]

This computes the square root of each element of the table using listability.

Sqrt[tab]

This computes timings for the two different ways of computing the square roots.  In each case the command is repeated 10000 times to get a better estimate.

{First[
        Timing[Do[
            Table[Sqrt[tab[[i]]],  {i,  
                Length[tab]}],  {10000}]]], 
First[
        Timing[Do[Sqrt[tab],  {10000}]]]}

Matrices

A List of vectors (Lists) of equal length represents a matrix.

Here is a list of three vectors.

A  =  {{1. ,  2,  3.   +  4.  I},  {1. ,  Sin[2],  3}, {7, 8, 
          9}}

It is a matrix.

MatrixQ[A]

You can view matrices in a grid form using MatrixForm

This shows the matrix with rows and columns lined up.

MatrixForm[A]

Matrices can also be generated by the Table command

This generates a matrix.

mat  =  Table[i^j, {i, 4}, {j, 4}]

Matrices can also be imported.

This loads a package of utilities for working with URLs.

Needs["Utilities`URLTools`"]

This creates a temporary file which holds the file referred to by the URL in uncompressed form.

tempfile  =  
      FetchURL["http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/\
bcsstruc1/bcsstk12.rsa.gz"]

This imports the file into Mathematica.

bc  =  Import[tempfile,  "HarwellBoeing"]

This makes a plot of the matrix structure.

ArrayPlot[bc@"ToPatternArray"[]];

Both Part and Take work for matrices (and more deeply nested objects)

This gets (Subscript[A, 21])

A[[2, 1]]

This gets  the first column of mat.

mat[[All, 1]]

This gets  lower left 2x2 submatrix from mat.

Take[mat,  -2, 2]

Listable commands work down to any level of nesting in lists.  This means that arithmetic operations on matrices work elementwise.

This divides each row of the matrix by its first element.

mat/{1. , 2. , 3. , 4. }

Linear Algebra

In Mathematica 5, the implementation of much of the numerical linear algebra was modified to use LAPACK (Linear Algebra PACKage), which is based on a small set of functions doing important but basic operations called the BLAS (Basic Linear Algebra Subroutines).  Besides robustness of the algorithms in LAPACK, the big advantage of this approach is that it is possible to get significant speed increases by using specially tuned BLAS libraries.

One of the basic BLAS operations is matrix multiplication.  The graph below shows the speedup one can get with a tuned BLAS library using multithreading versus simple compiled code.

There are many functions for doing different tasks in numerical linear algebra in Mathematica.  Consistent with the theme I mentioned at the beginning of this presentation, they are generally task based and automatically chose appropriate algorithms based on matrix type.

Because of their importance, I will demonstrate functions for solving systems of linear equations and for computing eigenvalues and eigenvectors.  For further reading I suggest Section 3.7 in The Mathematica Book and the extensive Advanced Documentation section on linear algebra.

The direct way for solving linear systems of equations from a matrix representation is LinearSolve.

This solves the system mat.x = b for b = {4,3,2,1} using exact rational arithmetic.

LinearSolve[mat,  {4, 3, 2, 1}]

When the matrix and other inputs are exact, exact arithmetic will be used and there will be no error.  Approximate arithmetic is much faster, but there can be error depending on matrix conditioning.  Matrices and other inputs will be coerced to the least common precision.

This solves the system mat.x = b for b = {4,3,2,1} using machine number arithmetic.

LinearSolve[mat,  N[{4, 3, 2, 1}]]

If you are going to solve a system repeatedly for different right hand sides, you can use LinearSolve with one argument, which typically computes a decomposition (usually LU) that can be used to solve for a given vector b quickly.

This sets up an LU decomposition for the matrix mat in a form which acts like an effective inverse.

lsf  =  LinearSolve[N[mat]]

This solves for a particular vector b = {1,2,4,8}

lsf[{1, 2, 4, 8}]

Just as for solutions of linear equations, the type of arithmetic used is based on the precision of the inputs.

This gives an exact expression for the eigenvalues.

Eigenvalues[mat]

This uses machine arithmetic to compute the eigenvalues.

Eigenvalues[N[mat]]

SparseArrays

When I imported the matrix from the Matrix Market web site previously, it was read in as a SparseArray.

bc

SparseArrays can be converted to dense matrices using Normal as long as sufficient memory is available.  The SparseArray explicitly stores only the values and positions in the matrix which are nonzero, so in this example it takes much less memory:

This compares the sizes of the two representations.

{ByteCount[bc],  ByteCount[Normal[bc]]}

SparseArray is also a convenient way to generate matrices by pattern.

This generates a differentiation matrix.

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

This shows it with rows and columns lined up.

MatrixForm[S]

Arithmetic for SparseArrays works by listability, so, just as for dense matrices, operations will be done elementwise.

This applies the differentiation matrix on a discretization of x (1 - x) on the unit interval.

h  =  1/11. ;
((S/(h^2))) . 
      Table[x (1  -  x),  {x,  h,  1  -  h,  h}]

Many operations that work for vectors and matrices work similarly for SparseArrays.

The most important tasks for numerical linear algebra work using sparse algorithms for SparseArray input.

This generates a convenient test vector.

b  = bc .  Range[Length[bc]];

This solves the system bc.x = b for x and shows the timing and part of the solution.

Short[Timing[LinearSolve[bc,  b]],  20]

There are also iterative algorithms, which can be chosen using the Method option.

This solves the system bc.x = b for x using a conjugate gradient algorithm and shows the timing and part of the solution.

Short[
      Timing[LinearSolve[bc,  b,  
          Method -> {Krylov,  Method -> ConjugateGradient,  
              Preconditioner -> "ILU"}]],  20]

For eigendecompositions and singular values, it is possible to specify that just a few of the characteristics be found, in which case a sparse algorithm can be used.

This finds the largest 2 eigenvalues of the matrix bc.

Eigenvalues[bc,  2]

Numerical Computations with Functions and Equations

In my view, the heart of numerical analysis is using numerical approximations to work with functions.

Mathematica has several numerical commands that work to approximate properties of functions such as FindMinimum, FindMaximum, NMinimize, NMaximize, NIntegrate, NSum, and NProduct.  There are also commands that work to numerically solve equations such as FindRoot, NSolve, and NDSolve.  Some commands such as FindFit and NonlinearRegress work with both data and functions.

There is not time to give you a reasonable demonstration of all of these, so I will demonstrate FindMinimum, which does local unconstrained optimization and NDSolve which solves differential equations with an emphasis on principles which are common to most of the numerical commands listed above.  The most important of these are:

Local Optimization

The command for finding local minima of functions is FindMinimum.  To demonstrate how it works, I will use a relatively simple function in 2-dimensions.

This loads a collection of standard test problems for unconstrained local optimization.

Needs["Optimization`UnconstrainedProblems`"]

The following minimization problem is the so-called "Rosenbrock"  problem, commonly used to check on local optimization methods in low dimensions.  It difficulty lies in a curved valley: simple descent methods do quite poorly on this problem.

p  =  GetFindMinimumProblem[Rosenbrock]

This is a good problem to compare algorithms since the 2-dimensions can easily be visualized.

The function FindMinimumPlot included in the UnconstrainedProblems package shows the steps and evaluations taken by FindMinimum to reach a solution.  The data for the plot is collected using the EvaluationMonitor and StepMonitor options.

This shows the steps and evaluations used by the method Mathematica chooses automatically.

fmp[1]  =  
      FindMinimumPlot[(((1 - (Subscript[X, 1])))^2) + 
          100 (((-(Subscript[X, 1]^2) + (Subscript[X, 2])))^2), {{(Subscript[X, 1]), -1.2`}, {(Subscript[X, 2]), 
            1.`}}]

This shows the steps and evaluations used by the Newton's method .

fmp[2]  =  
      FindMinimumPlot[(((1 - (Subscript[X, 1])))^2) + 
          100 (((-(Subscript[X, 1]^2) + (Subscript[X, 2])))^2), {{(Subscript[X, 1]), -1.2`}, {(Subscript[X, 2]), 
            1.`}},  Method -> Newton]

This shows the steps and evaluations used by the BFGS quasi-Newton method .

fmp[3]  =  
      FindMinimumPlot[(((1 - (Subscript[X, 1])))^2) + 
          100 (((-(Subscript[X, 1]^2) + (Subscript[X, 2])))^2), {{(Subscript[X, 1]), -1.2`}, {(Subscript[X, 2]), 
            1.`}},  Method -> QuasiNewton]

All of the methods above use derivatives which are computed symbolically by Mathematica.  When Mathematica cannot compute derivatives symbolically, such as when given a black-box function, finite differences are used instead.

This shows the steps and evaluations used by the BFGS quasi-Newton method with derivatives computed by finite differences.

fmp[4]  =  
      FindMinimumPlot[(((1 - (Subscript[X, 1])))^2) + 
          100 (((-(Subscript[X, 1]^2) + (Subscript[X, 2])))^2), {{(Subscript[X, 1]), -1.2`}, {(Subscript[X, 2]), 
            1.`}},  Method -> QuasiNewton,  
        Gradient -> "FiniteDifference"]

There is also a derivative-free method build into FindMinimum, which is based on the Principal Axis method of Brent.  It needs two starting conditions to compute secant derivatives to get started.  When you give two starting conditions to FindMinimum, it uses this algorithm.

This shows the steps and evaluations used by the BFGS quasi-Newton method with derivatives computed by finite differences.

fmp[5]  =  
      FindMinimumPlot[(((1 - (Subscript[X, 1])))^2) + 
          100 (((-(Subscript[X, 1]^2) + (Subscript[X, 2])))^2), {{(Subscript[X, 1]), -1.2`,  -1. }, {(Subscript[X, 2]), 1.`,   .5}}]

Below is a graph and table which summarizes.

({{"", "LevenbergMarquardt", "Newton", "BFGS", "FD BFGS", "PrincipalAxis"}, {"", 0.`, 4.9696177196758636`*^-18, 3.009122763698238`*^-21, 2.002901521029204`*^-11, 1.9162112723944665`*^-11}, {"", (Subscript[X, 1]) -> 1.`, (Subscript[X, 1]) -> 
                0.9999999977999193`, (Subscript[X, 1]) -> 
                1.0000000000278046`, (Subscript[X, 1]) -> 
                0.9999955249004842`, (Subscript[X, 1]) -> 1.000002403410144`}, {"", (Subscript[X, 2]) -> 1.`, (Subscript[X, 2]) -> 
                0.9999999955638855`, (Subscript[X, 2]) -> 
                1.0000000000508806`, (Subscript[X, 2]) -> 
                0.9999910448214613`, (Subscript[X, 2]) -> 1.0000044409608926`}, {"", "Steps" -> 15, "Steps" -> 
                22, "Steps" -> 42, "Steps" -> 
                41, "Steps" -> 15}, {"", "Residual" -> 21, "Function" -> 
                29, "Function" -> 51, "Function" -> 
                216, "Function" -> 596}, {"", "Jacobian" -> 16, "Gradient" -> 
                29, "Gradient" -> 51, "Gradient" -> 
                54, ""}, {"", "", "Hessian" -> 
                23, "", "", ""}})

Another interesting feature of the local optimization functions is the adaptive precision control which is above machine number precision.  You can specify the maximum precision that will be used by using the WorkingPrecision option.  Default tolerances are chosen to be about half of this precision.

This uses Newton's method to find the minimum of the Rosenbrock function using adaptive precision up to 100.  The distance from the solution and current precision are printed at each step.

FindMinimum[(((1 - (Subscript[X, 1])))^2) + 
        100 (((-(Subscript[X, 1]^2) + (Subscript[X, 2])))^2), {{(Subscript[X, 1]), -1.2}, {(Subscript[X, 2]), 
          1. }},  
Method -> 
        Newton, 
WorkingPrecision -> 
        1000, 
StepMonitor :> 
        Print[{N[Norm[{(Subscript[X, 1]),  (Subscript[X, 2])}  -  {1, 1}]], 
            Precision[{(Subscript[X, 1]),  (Subscript[X, 2])}]}]]

Differential Equations

The command for solving differential equations is NDSolve.  In addition to

It also features

One of the most useful features of NDSolve for scientists is its ability to solve equations of different type using appropriate algorithms.  A good example of this is a problem in chemical kinetics due to Robertson.

kinetics  =  {(Subscript[y, 1])'[
              t]  == -(1/25) (Subscript[y, 1])[t]  +  (10^4) 
                (Subscript[y, 2])[t] (Subscript[y, 3])[t],  (Subscript[y, 2])'[t]  ==  
            (1/25) (Subscript[y, 1])[t]  -  (10^4) 
                (Subscript[y, 2])[t] (Subscript[y, 3])[t] -   3 (10^7) 
                (Subscript[y, 2])[t]^2 };
balance  =  
        (Subscript[y, 1])[t]  +  (Subscript[y, 2])[t]  +  (Subscript[y, 3])[t]  ==  
          1;
start  =  {(Subscript[y, 1])[0]   ==  1,  (Subscript[y, 2])[0]  ==  0,  
          (Subscript[y, 3])[0]  ==  0};

This system is often solved by differentiating the balance equation to obtain a system of three ODEs.

This solves the Robertson system as a system of three ODEs.

sol  =  
      NDSolve[{kinetics,  D[balance,  t],  start},  {(Subscript[y, 1]),  (Subscript[y, 2]),  
          (Subscript[y, 3])},  {t, 0, (10^6)}]

The solution is returned is effectively a set of three functions.  InterpolatingFunctions interpolate between the data stored by the method used to compute the solution to get intermediate values when requested.

This loads two packages which will be helpful in making an informative plot of the solutions.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["Graphics`Graphics`"];

This makes a plot of the solutions with t on a logarithmic scale. (Subscript[y, 1])is shown in red, (Subscript[y, 2]) is scaled by a factor of (10^4)and is shown in green, (Subscript[y, 3]) is shown in blue.  The points shown are at the steps that the method used to compute the solution.

times  =  
        Drop[First[InterpolatingFunctionCoordinates[(Subscript[y, 1]) /.  First[sol]]], 
          1];
Show[
        Block[{$DisplayFunction  =  Identity}, {LogLinearListPlot[
              Transpose[{times,  (Subscript[y, 1])[times]  /.  First[sol]}], 
              PlotStyle -> Red], 
LogLinearListPlot[
              Transpose[{times,  
                  (10^(4 )) (Subscript[y, 2])[times]  /.  First[sol]}], 
              PlotStyle -> Green], 
LogLinearListPlot[
              Transpose[{times,  (Subscript[y, 3])[times]  /.  First[sol]}], 
              PlotStyle -> Blue]}],  
        PlotRange -> {{-6, 6}, {0, 1}}];

You can see from the plot that (Subscript[y, 2]) has varied on a much quicker time scale than (Subscript[y, 1])or (Subscript[y, 3]).  A problem which has different time scales like this that can affect solver stability is known as stiff.  NDSolve automatically switches to a stiff solver when appropriate.

Note that if you try to solve this system with an explicit method, it will get bogged down quickly.

sol  =  
      NDSolve[{kinetics,  D[balance,  t],  start},  {(Subscript[y, 1]),  (Subscript[y, 2]),  
          (Subscript[y, 3])},  {t, 0, (10^6)},  
        Method -> ExplicitRungeKutta]

Another way to solve the problem is to leave the balance as an algebraic condition and solve the system of differential-algebraic equations (DAEs).  When you do this, NDSolve automatically recognizes the equations as DAEs and uses a method designed for DAEs.

This solves the Robertson system as a system of DAEs.

daesol  =  
      NDSolve[{kinetics,  balance,  start},  {(Subscript[y, 1]),  (Subscript[y, 2]),  (Subscript[y, 3])},  {t, 
          0, (10^6)}]

daetimes  =  
        Drop[First[
            InterpolatingFunctionCoordinates[(Subscript[y, 1]) /.  First[daesol]]], 
          1];
berr  =  Apply[Subtract,  balance];
Show[
        Block[{$DisplayFunction  =  Identity}, {LogLogListPlot[
              Transpose[{times,  
                  Abs[(berr  /.  First[sol])  /.  t -> times]}], 
              PlotStyle -> Red], 
LogLogListPlot[
              Transpose[{daetimes,  
                  Abs[(berr  /.  First[daesol])  /.  
                      t -> daetimes]}], PlotStyle -> Green]}],  
        PlotRange -> {{-6, 6}, {-16, -14}}];

Something new in version 5.1 is an interactive interface for planar systems called EquationTrekker.

This loads the EquationTrekker package.

Needs["NumericalMath`EquationTrekker`"]

This opens a window which allows interactive selection of initial conditions and parameters for a harmonic oscillator with nonlinear driving.

EquationTrekker[x''[t]  +  
          x[t]  ==  α Cos[ω x[t]], 
      x, {t, -10, 10},  
      PlotRange -> {{-10, 10}, {-10, 10}}, 
      TrekParameters -> {α -> 3,  ω -> 
            1}]

This allows interactive selection of starting condition and parameters for the Poincare map of the Duffing oscillator.

\(EquationTrekker[{\ \(x''\)[t]\  + \ γ\ \(x'\)[\ t]\  - \ 
            x[t]\  + \ x[t]^3\  == \ ϵ\ Cos[ω\ t]}, \ 
      x, \ {t, 1000, \ 
        10000}, 
PlotRange -> {{\(-2\), 2}, {\(-2\), 
            2}}, 
TrekParameters -> {γ ->  .15, \ ϵ ->  .3, \ ω -> 1. }, \ 
      TrekGenerator -> {PoincareSection, \ "SectionCondition" -> Mod[ω\ t, \ 2\ π], \ "SectionVariables" -> {x, \ 
              x'}, MaxSteps -> ∞}]\)

NDSolve is also able to solve a limited set of PDEs: 1 + n dimensional evolution equations which have initial conditions in the time like direction.

Differential Equations

For nearly all of the numerical commands that work with functions, error tolerances are quite easy to specify using the AccuracyGoal and PrecisionGoal options.  The exact interpretation depends on the command being used, but generally, PrecisionGoal->(Subscript[G, p]) and AccuracyGoal->(Subscript[G, a]) limits ‖(Subscript[err, X])‖ < (10^(-(Subscript[G, a]))) + ‖X‖(10^(-(Subscript[G, p]))), where (Subscript[err, X ])is the estimated error in the approximate solution X.  (10^(-(Subscript[G, a]))) is the absolute error tolerance and (10^(-(Subscript[G, p]))) is the relative error tolerance.

Often, if you want very strict tolerances, you may need to increase the WorkingPrecision to achieve them.  For many operations with functions if you have WorkingPrecision->p, it may be unrealistic to, in general, expect to achieve precision goals (Subscript[G, p])> p/2.  Many of the Mathematica commands have defaults (Subscript[G, p])  =  p/2 and (Subscript[G, a])  =  p/2 .

As a final example, I'll show you some solutions with NDSolve using different tolerances as an example.

This solves the pendulum equation with default tolerances.

Timing[
      sd  =  First[
          NDSolve[{x''[t]  + Sin[x[t]]  ==  0,  
              x[0]  ==  3.14159,  x'[0]  ==  0}, 
            x, {t, 0, 1000}]]]

This shows a phase-plane plot of the solution near the saddle point on the separatrix.

ParametricPlot[{x[t],  x'[t]}  /.  sd, {t, 0, 1000},  
      PlotRange -> {{3.14,  3.142}, {- .002,  .002}}, 
      PlotPoints -> 10000]

Clearly from the plot, there is some noticeable error.

This solves the pendulum equation with 12 digit tolerances.

Timing[
      sd  =  First[
          NDSolve[{x''[t]  + Sin[x[t]]  ==  0,  
              x[0]  ==  3.14159,  x'[0]  ==  0}, 
            x, {t, 0, 1000},  AccuracyGoal -> 12,  
            PrecisionGoal -> 12,  
            MaxSteps -> Infinity]]]

This shows a phase-plane plot of the solution near the saddle point on the separatrix.

ParametricPlot[{x[t],  x'[t]}  /.  sd, {t, 0, 1000},  
        PlotRange -> {{3.1415,  3.1416}, {- .0002,  .0002}}, 
        PlotPoints -> 10000];

Much less error, but it is still there.  To get much better tolerances, you need to go to higher precision.

This solves the pendulum equation with 16 digit tolerances.  To achieve these tolerances it is necessary to set the working precision to a sufficiently high value.

Timing[
      sd  =  First[
          NDSolve[{x''[t]  + Sin[x[t]]  ==  0,  
              x[0]  ==  
                SetPrecision[3.14159,  32],  x'[0]  ==  0}, 
            x, {t, 0, 1000},  AccuracyGoal -> 16,  
            PrecisionGoal -> 16,  MaxSteps -> Infinity,  
            WorkingPrecision -> 32,  
            Method -> Extrapolation]]]

This shows a phase-plane plot of the solution near the saddle point on the separatrix.

ParametricPlot[{x[t],  x'[t]}  /.  sd, {t, 0, 1000},  
        PlotRange -> {{3.1415899,  
              3.1415901}, {- .000001,  .000001}}, 
        PlotPoints -> 10000];

In this talk, I have taken you through a whirlwind introduction to numerical computing with Mathematica.

Support for a variety of different types of numbers with consistent rules for arithmetic and numerical coercion form the foundation for numerical computing in Mathematica.

Vectors, matrices, and commands supporting them form the next level of the base for numerical computing in Mathematica.  Strong support for numerical linear algebra on both sparse and dense matrices is important for higher level commands.

The symbolic and strong numerical base make it possible for Mathematica to have a powerful set of commands for doing numerical operations on functions.