Numerical Computation in Mathematica

Mark Sofroniou and Rob Knapp/Wolfram Research

Mark has shown some of the nuts and bolts of the system along with NDSolve.  I’ll focus on numerical procedures that work with functions.   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.

General Design Principles for 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, NIntegrate instead of Newton, Simpson).

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.

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

The default is generally to use machine precision.  In this case, functions are compiled to improve evaluation speed.

Error is controlled by options AccuracyGoal and PrecisionGoal.

In general, the algorithm works until |e|≤ max((10^(-(Subscript[g, a]))), r (10^(-(Subscript[g, p])))) where e is the estimated error and r is the approximate result.

The setting for WorkingPrecision needs to be high enough so that these goals can be satisfied

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

Numerical Integration

This numerically integrates a highly oscillatory function

NIntegrate[(1/x)Cos[Log[x]/x],{x,0,1}]

Plot[(1/x)Cos[Log[x]/x],{x,0,1}, PlotRange->20]

You can see where the function was evaluated using EvaluationMonitor:

{r, {evalx}} = Reap[NIntegrate[(1/x)Cos[Log[x]/x], {x,0,1}, EvaluationMonitor:>Sow[x]]]; Column[{ListPlot[Transpose[{evalx, Range[Length[evalx]]}], ImageSize->Large],ListPlot[Transpose[{x,  (1/x)Cos[Log[x]/x]} /. x->evalx],Filling->0, ImageSize->Large]}]

Surely there is some magic going on.  One way to check is to repeat, but setting up the function as a “black box”

SetAttributes[bb, Listable];bb[x_?NumberQ] := (1/x)Cos[Log[x]/x];

{r, {evalx}} = Reap[NIntegrate[bb[x], {x,0,1}, EvaluationMonitor:>Sow[x]]]; Column[{ListPlot[Transpose[{evalx, Range[Length[evalx]]}], ImageSize->Large],ListPlot[Transpose[{x,  (1/x)Cos[Log[x]/x]} /. x->evalx],Filling->0, ImageSize->Large]}]

What is happening is that NIntegrate is able to make transformations based on the symbolic form of the function so that it can use a Levin rule method.

PrecisionGoal

Consider a problem with an exact solution:

exact = Integrate[1/Sqrt[x+y],{x,0,1},{y,0,x}]

AbsoluteTiming[default = NIntegrate[1/Sqrt[x+y],{x,0,1},{y,0,x}]]

exact - default

Increase the precision goal.

AbsoluteTiming[better = NIntegrate[1/Sqrt[x+y],{x,0,1},{y,0,x}, PrecisionGoal->12]]

With PrecisionGoal->12, the error should be below

better (10^(-12))

better - exact

AbsoluteTiming[attempt = NIntegrate[1/Sqrt[x+y],{x,0,1},{y,0,x}, PrecisionGoal->18]]

attempt - exact

AbsoluteTiming[best = NIntegrate[1/Sqrt[x+y],{x,0,1},{y,0,x}, PrecisionGoal->18, WorkingPrecision->40]]

exact - best

Local Optimization

Local optimization functions FindMinimum, etc start at a search point and try to find a local optimal point.

m1 = FindMinimum[Cos[x] - Sin[x],x]

The starting point can affect the result

m2 = FindMinimum[Cos[x] - Sin[x],{x,-.56}]

{m1[[1]] - m2[[1]], (x /. m1[[2]]) - (x /. m2[[2]])}

Lets see why the message was issued.

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

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

One solution is to use higher precision

FindMinimum[Cos[x] - Sin[x],{x,.56}, AccuracyGoal->50, PrecisionGoal->50, WorkingPrecision->150]

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.

s[x_, y_] = Cos[2+9 (((x-y))^2)]+Cos[1+7 (y^2)]+Sin[5 (x^2)]+Sin[3+4 (((x+y))^2)];
ps = 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.

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

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

NMinimize will work harder to attempt to find a global minimum.  This directs it to use the random search method.

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[ps, Graphics3D[{Red, PointSize[0.05],Point[{x, y, lm} /. lmxy]}],Graphics3D[{Green, PointSize[0.05],Point[{x, y, gm} /. gmxy]}]]

Final example: Fitting to an integrated value.

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

Here is some data:

ListPlot[data]

For a model function use

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

Note how PatternTest was used to prevent symbolic evaluation.  In this case, symbolic evaluation would only result in messages being issued, but there are many instances where symbolic evaluation effectively changes the function,

Now use FindFit to find the best fit.  The PrecisionGoal and AccuracyGoal were relaxed because the model only has limited accuracy,

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

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

That worked, but it was not very efficient.  A different way to do this is to use NDSolve for the model:

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

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

This is still much more inefficient than it might be.  The problem is that NDSolve is being evaluated for each separate value of x,  Since NDSolve can compute all the values on [0, 4] for a given a, b it is better to compute the solution just once for any value {a, b}.  This can be done by adding a rule to cache the values in the definition

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

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

You may want to Clear the definition since the saved values may take some memory

DownValues[model1]

The end

Mark and I have given you a whirlwind tour of some of the numerical computation features in Mathematica, but have left so many things untouched.  I highly recommend looking at the documentation for examples and information.  Some useful more advanced tutorials in the collection are:

Advanced Numerical Differential Equation Solving in Mathematica

Advanced Numerical Integration in Mathematica

Unconstrained Optimization

Constrained Optimization