Mathematica Numerics

Mathematica has many powerful built-in numerical algorithms.  This talk will provide an introduction to the effective use of algorithms for numerically computing, roots, minima, integrals, sums, solutions of differential equations, and matrix quantities.

Initialization

Needs["Graphics`Colors`"];

Needs["Graphics`Graphics`"];

Off[General::spell];

Off[ParametricPlot::ppcom];

eps = 10.^(-$MachinePrecision)

A Tour of Numerical Functions in Mathematica

If all problems could be solved analytically, there would be no need for numerical approximations.  However, there are many problems at the cutting edge of research and in applications for which there is no analytic solution.  The numerical functions in Mathematica offer a variety of sophisticated numerical algorithms that allow you to continue on when symbolic calculations become difficult or impossible.

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

Root finding.

FindRoot[]

Of course, the square root function is built into Mathematica, so it is not necessary to use FindRoot[] to compute the square root of 2, but it provides a good example illustrating some aspects of the function.

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

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

Now suppose we wanted to compare the rate of convergence of Newton's method with that of the secant method.  It would be nice to see the sequence of points used.  This can be done relatively easily by

FindRoot[Print[InputForm[x]];x^2 == 2,{x,1}]

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

How does this work?  Usually with Mathematica commands, the arguments are evaluated, as in

D[Print[InputForm[x]];x^2-2,x]

However some commands, including the command FindRoot[] have the attirbute HoldAll.  This tells Mathematica to leave the arguments unevaluted.  When FindRoot forms the function which is to be evaluated at the values of x, it keeps the entire function, thus each time a new value of x is tried, the value is printed.

Many of the numerical functions in Mathematica have the HoldAll attribute.  If you are not sure, you can always check using the Attributes[] command:

Attributes[FindRoot]

It is interesting to note that Newton's method required 1 symbolic evaluation (to find the derivative), but it converged much more quickly to the result.

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.

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

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},
	Contours->{1,2},ContourShading->False];

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

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

A plot also verifies that we have the right solution:

Plot[{t,b^t /. %},{t,0,4},PlotStyle->{Red,Black}];

NSolve

FindRoot has the advantage that it can solve transcendental equations.  However, it has the serious drawback that it only finds a single root.  If you want to find all the roots of an algebraic equations, it is better to use NSolve.

NSolve[x^2 == 2,x]

Keep in mind, though that NSolve will not work for many equations, such as the problem we successfully solved with FindRoot above:

NSolve[{b^x == x,D[b^x,x] == 1},{x,b}]

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[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[Sin[x]],{x,1,10000}]]

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

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

Wow!  That NSum was way faster!  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 which 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}]

Integration

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 reagions separately by including the points where there might be discontinuities.

Plot[f[x],{x,-3,6}];

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

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 most any Mathematica function which evaluates to a number.

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 which vary on extremely different time scales.  For example,

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

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

Linear Algebra

Mathematica has symbolic as well as numerical linear algebra capabilities.  Generally for large matrices, the symbolic calculations are impracticaland numerical computations become important.  Here is an example of some of the operations you can do with a matrix.

First construct a square matrix of random real values:

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

Since the matrix was constructed randomly, you don't know apriori 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.

Minv = Inverse[M]

We can check on the inverse by taking the matix product, which is denoted by "." in Mathematica.

M.Minv

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

Chop[%] // MatrixForm

Define a random vector b

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

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

x1 = Minv.b

A better 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:

{evals,evecs} = Eigensystem[M]

Chop[M-Transpose[Inverse[evecs].
	 			 DiagonalMatrix[evals].
		  		 evecs]]

Transpose[Inverse[evecs].
		  DiagonalMatrix[Exp[evals]].
		  evecs]

MatrixExp[M]

Combinations

One of the powers of Mathematica is that you can easily put together combinations of commands.  For example, the following program can be used to solve a boundary value problem in ODE's:

Clear[f];
f[lambda_?NumberQ] := NDSolve[{
	y''[x] + lambda^2 y[x] == 0,
	y[0] == 0,y'[0] == 1},
	y,{x,0,1}][[1,1,2]][1]

When you do this, you have to be thoughful of what is being done.  For example,

FindRoot[f[lambda],{lambda,1}]

Does not work because Mathematica cannot find the derivative of f[lambda] with respect to lambda.  However, by specifying two starting points,

FindRoot[f[lambda],{lambda,1,2}]

You can find other eigenvalues by using different starting values

FindRoot[f[lambda],{lambda,4,10}]

Here is a way to get all of the eigenvalues between 0 and 20

plot = Plot[f[lambda],{lambda,0,20}];

A part of the plot is a list containing the points.  This list can be isolated from the rest of the plots.

Position[plot,Line[___]]

list = plot[[1,1,1,1]];

Now you can make a list of starting points.  Save the lambda coordinates of the points where the sign of f[lambda] changes.

starts = {};
Do[If[Sign[list[[j,2]]] != Sign[list[[j+1,2]]],
	starts = Join[starts,{{list[[j,1]],list[[j+1,1]]}}]],
	{j,1,Length[list]-1}];
starts

starts = {{2.916666666666667, 3.333333333333333}, 
 {6.25, 6.666666666666666}, 
 {9.16666666666667, 9.58333333333333}, 
 {12.5, 12.91666666666667}, 
 {15., 15.83333333333334}, 
 {18.75, 19.16666666666667}}

Now using Map[], let FindRoot get the eigenvalues.

eigenvals = (lambda /. 
FindRoot[f[lambda],{lambda,#}])& /@ starts

Now you can use the NDSolve command to obtain the eigenfunctions:

eigenfuncs = (NDSolve[{
	y''[x] + #^2 y[x] == 0,
	y[0] == 0,y'[0] == 1},
	y[x],{x,0,1}][[1,1,2]])& /@ eigenvals;

Plot[Evaluate[eigenfuncs],{x,0,1},PlotRange->All]

Of course, the next step might be to make them an orthonormal set with respect to the inner product!

ip[f1,f2,x] := Integrate[f1*f2,{x,0,1}]

Accuracy and Precision in Arithmetic

In Mathematica, there are essentially two different types of numerical quantities, exact and approximate.  Naturally, exact quantities have the advantage that computations done with them have no numerical error.  At the same time,  calculations done with exact numbers can be slow and in many cases can only be done symbolically, and so, often, the results are too large to be understandable.  .  Generally, if you need to have a result, approximate numbers will give you a more compact result much more quickly.  However, there may be error in the approximate results.

Example: Coefficients with respect to a basis

For example, consider the problem of representing vectors in the plane with repect to a basis.  If the basis vectors are nearly orthogonal, there is no problem.  However, if the basis vectors are in nearly the same direction, numerical errors may be introduced.

Let the matrix of basis vectors be

B = {{12/13,12/13},{9999999999,10000000001}};

The basis vectors are very nearly in the same direction.  As a consequence, large multiples of the vectors will have to be added to get vectors in other directions.

Compute the coefficents of the vector {E,Pi} with repect to the basis B.

coefs = LinearSolve[B,{E,Pi}]

This computation has been done in exact arithmetic, using integers, rationals, and symbolic representations of the numbers E and Pi.

Since the computation was done exactly, constructing the vector from the basis coefficients should also work out exactly.  It does, once it is simplified

Simplify[B.coefs]

With more complicated expressions, the simplification may not put the expression into a form simple enough to be understood.

Now try to construct the vector using numerical approximations of the coefficients

B.N[coefs]

Even though N[coefs] gives an approximation of the coefficients with  machine precision (16 digits of precision on this machine),  the result clearly has no digits of precision.  What has happened is that even though the errors in machine numbers are small, those can be magnified by multiplication.  When the numbers with magnified error are later subtracted, the error may be left over.

Bignums: Precision and Accuracy

In many cases, however, it is not possible to do the compuation with exact numbers, either because the results do not simplify or because the computation would take too much time.  In addition to machine numbers, Mathematica provides another inexact approximation to numbers, called bignums.  Bignums are a compromise between the speed of machine numbers and the accuracy of exact computations.  In addition to doing basic arithmetic operations to an arbitrary specified precision, bignums keep track of what the error in the result may be and adjusts the precision accordingly.

Now try the computation with a bignum approximation to the coefficients with 25 digits of precision.

B.N[coefs,25]

Mathematica shows the number of digits it believes are correct based on the arithmetic.

Precision[]

The command Precision[] allows you to find out the number of digits of precision are in an arithmetic result.

Precision /@ %

Precision /@ (B.N[coefs])

Does this mean that the machine number calculation got a better result for the first value?  The answer is no.  Because machine arithmetic is done with the hardware, Mathematica cannot keep track of the error possible in the results.  For this reason, for all machine numbers except machine zero (0.) the result of Precision[] will be equal to $MachinePrecision.

While the bignum arithmetic is much slower than machine arithmetic, the way it keeps track of its error can be a tremendous advantage--it gives you a worst case estimate of the arithmetic error in a result.

Precision gives a relative measure of how good your numerical approximation is.  The measure of precision is essentially the numbers of digits you know correctly whatever the size of the numbers.  However,  the notion of precision does not always contain the information you wnat about error.  For example,

B.N[coefs,25] - N[{E,Pi},25]

has precision

Precision /@ %

Accuracy[]

It may be the case that no digits are really known, but that does not mean that your answer is totally meaningless.  What you do know in this case even if you don't even know the first digit is the size of the result, or in another sense the  absolute error in the arithmetic result.  Mathematica provides the function Accuracy[] to measure just this quantity.

Accuracy /@ (B.N[coefs,25] - N[{E,Pi},25])

Which gives you the information that even though the result may have zero precision, you know both terms to within an absolute error of 10^-4, which makes sense since this is how many digits should have cancelled based on the precision of your results.

Scale

Precision is the number of digits in the number itself known to be be correct.  Accuracy is the number of digits which are known to the right of the decimal point--the absolute error in the result.  Let the scale of a number be the logarithm to base 10 of the numbers absolute value, or the number of digits to the left of the decimal place.  Then, it should be that 

	Precision - Accuracy = Scale
	
for any number.  Some examples will help illustrate this relationship and improve your understanding of the distiction between accuracy and precision.

Clear[f];
f[x_] := {Precision[x],Accuracy[x],
	Log[10,Abs[x]],Precision[x]-Accuracy[x]}

f[N[1,20]]

f[N[10000000000000000000000000,20]]

f[N[1/10000,20]]

f[1+N[10^-100,20]]

How did you get precision 120 here?  When adding to exact numbers, the precision will be increased up to a certain point to obtain a good result.

f[1+N[10^-100000,20]]

The maximum precision Mathematica will go to is determined by the variable

$MaxPrecision

SetPrecision[] and SetAccuracy[]

It may be the case that the precision or accuracy that Mathematica returns with a result is not appropriate for further calculations that you need to do.  In this case, there are two commands to change the precision and accuracy of numbers, SetPrecision and SetAccuracy.

SetPrecision[B.N[coefs,25],25]

When Mathematica computes the precision and accuracy lost with arithmetic operations, it uses a "worst case" estimate.  In other words, the precision shown is the number of digits that are certain to be accurate.  It is often the case that there are actually more digits which do happen to be accurate.  In the vector example:

% - N[{E,Pi},25]

there were quite a few more correct digits than the result of Precision[] indicated.

You can use SetPrecision to set a number to any precision you like.  If you choose a higher value of precision than the number has digits, Mathematica will add incorrect digits to the number.  For example,

third = N[1/3,17]

SetAccuracy[third,50]

You may notice that there were more correct digits than 17.  This is because Mathematica keeps a few extra digits to ensure that the cacluations are accurate to the precision and accuracy stated.  If you use SetPrecision or SetAccuracy to lower the number of digits, additional digits will be truncated from the number and the result will be rounded appropriately.

Or if you can specifiy infinite accuracy, Mathematica returns an exact number (almost always a rational) which has a decimal expansion which matches all of the digits of accuracy in the inexact number it came form.

SetAccuracy[third,Infinity]

Not quite the fraction you had in mind.  Often, the function Rationalize does more what you would like for cases like this:

Rationalize[third]

Example: Least Squares Fit to data.

The example shown above with two basis vectors which have very nearly the same direction may seem contrived, but the situation actually comes up quite often in numerical calculations.  Below is an example which demonstrates this in a more realistic setting.

Suppose we have 16 pairs of data points,

data = Table[N[{x,Sin[2 Pi x]}],{x,10,11,1/15}]

Suppose that we would like to find an interpolating polynomial of degree 15 which passes through all 16 points.  One way to do this is to do a fit.  Since there is just enough data to determine each coefficient in the polynomial, the error at each data point will be zero and thus it produces an interpolating polynomial.  To use the Mathematica function Fit[], we need to define a set of basis functions.  One choice is

basis = x^Range[0,15]

poly1[x_] = Fit[data,basis,x]

It is worth checking to see if the polynomial actually passes through the points.

(poly1[#[[1]]] - #[[2]])& /@ data

Numerical roundoff error has crept in and made the result quite poor.  What has happened is that some of the basis functions are barely distinguishable on the interval [10,11].  For example, if we normalize x^14 and x^15 to have the same value 1 at x = 10, we can compare the difference in the shapes of the functions

Plot[{(x/10)^14,(x/10)^15},{x,10,11}];

The functions do not differ by much on this interval.

The problem is not with Fit[], but in the poor choice of basis functions.  Generally, you will be better off using a set of basis functions which will be significantly different on the interval in question..  For example the result can be dramatically improved by using powers of x-10 instead of powers of x.

shiftbasis = Table[(x-10)^j,{j,0,15}];

poly2[x_] = Fit[data,shiftbasis,x];

(poly2[#[[1]]] - #[[2]])& /@ data

An even better way to solve this problem is to use the function InterpolatingPolynomial[], which computed the interpolating polynomial in a form so that evaluation of it minimizes numerical error.

poly3[x_] = InterpolatingPolynomial[data,x]

(poly3[#[[1]]] - #[[2]])& /@ data

Notice that much of the improvement in accuracy comes for the way in which the polynomial is evaluated. If you expand it  first, then you get

poly4[x_] = Expand[poly3[x]];
(poly4[#[[1]]] - #[[2]])& /@ data

In summary, then there are two types of numbers, exact and inexact.   Further, there are two types of inexact numbers, machine numbers and bignums.  Machine numbers are fast, but do not contain any information about the precision or accuracy of arithmetic results.  Bignums can be given whatever precision you like and contain information which can tell you about the accuracy and precision of arithmetic results.

Accuracy and Precision in Numerical Calculations

The following numerical functions take the options AccuracyGoal and PrecisionGoal (with the exception of FindRoot).  The default values for these options is given in the following table. (Numbers in parentheses represent the actual value with the default value of the WorkingPrecision option, which uses machine numbers)

Function		AccuracyGoal	PrecisionGoal

FindMinimum		Automatic (6)		Automatic (6)
FindRoot		Automatic (6)		N/A
NIntegrate		Infinity			Automatic (6)
NDSolve		Automatic (6)		Automatic (6)
NProduct		Infinity			Automatic (6)
NSum			Infinity			Automatic (6)

In each case, these options affect the numerical accuracy and precision of the result.  Note that the word "affect" was used as opposed to "determine"--the examples below show that the accuracy and precision achieved depends very much on the problem being solved.  Including the word "Goal" in the name of the options is most appropriate.

The general model for error control

The general model used is

ME <= 10^(-AccuracyGoal) + MS 10^(-PrecisionGoal)

where ME is some measure (usually an estimate) of the error in the calculation, and MS is some measure of the size, or magnitude, of the solution.

In effect, the AccuracyGoal option represents the accuracy, or absolute, error tolerance in decimal digits, just as the function Accuracy[a] returns the number of digits to the right of the decimal point in the number a -- i.e. a measure of the absolute error in the number.

Similarily,  the PrecisionGoal option represents the relative error  tolerance in decimal digits, just as the function Precision[a] returns the total number of significant digits in the number a.  i.e. a measure of the relative error in the number.

It is extremely important to note that the in error model described above, the absolute tolerance from AccuracyGoal and the relative tolerance from PrecisionGoal are added together, providing a less stringent error tolerance than either one alone. 

 This is appropriate to have in a default error control model because it is not known a priori what the magnitude of the solution will be.  For solutions with MS >> 1, absolute error is overshadowed by the relative error.  Similarily, if MS is near 0, relative error may prove much too stringent to allow for a realistic solution.  
 
 If it is necessary to control one of absolute error or relative error, this is easily achived by setting PrecisionGoal or AccuracyGoal to Infinity, respectively.  With an infinite exponent, the corresponding term is zero, and the desired one remains.

AccuracyGoal: FindRoot

The function FindRoot has the option AccuracyGoal, but not the PrecisionGoal option.  This is due to the object of the function:  to find a root.  Near zero, relative error is ambiguous, so absolute error is the only appropriate measure.

FindRoot continues iterating until the absolute value of the objective function is less than or equal to the absolute tolerance specified by the value of the AccuracyGoal option.  This is not always what users of the function expect, as the following example shows.

mypi = x /. FindRoot[Cos[x]+1,{x,3}]

The default value of the AccuracyGoal option with machine precision calculations is 6.  The solution certainly doesn't seem to be within 10^-6 of Pi.

N[Pi-mypi]

In fact it is not even within 10^-3!  FindRoot has not, however, failed to do its task.  It is expected to find a value of the variable so that the absolute value of the objective function meets the tolerance specified by the AccuracyGoal option:

Abs[Cos[mypi]+1]

The problem here is that asymptotically, Cos[x]+1 looks very much like (x-Pi)^2, which has a quadratic tangent root -- the size of the objective function decreases quadratically for a linear change in the variable as it approaches the root.  A plot showing the linear approximations used by two succesive steps of Newton's method shows that each step makes it about half of the remaining distance to the root.

Show[Plot[Cos[x]+1,{x,3,3.25},DisplayFunction->Identity],
	Graphics[{Red,Line[{
		{3,Cos[3]+1},
		{3.070914844302653,0}}]}],
	Graphics[{Green,Line[{
		{3.070914844302653,Cos[3.070914844302653]+1},
		{3.106268467156334,0}}]}],
	DisplayFunction->$DisplayFunction];

It is clear from the plot that if the steps were doubled, the convergence would be much faster.  FindRoot allows you to add a factor to the size of each step using the option DampingFactor.  For problems like this where the root is quadratic, teh setting DampingFactor->2 will dramatically speed up convergence.

A good reason to use the size of the objective function as measure of the error is that it is one of the few situations in numerical calculations where the error is known exactly!  Because the error is known exactly, the algorithm can be made more robust.

It is occasionally asked what one should do if you need to know that the root is within some tolerance as opposed to the objective function.  The following is a crude technique for achieving this goal.

A function to attach the error estimate to the approximate distance to the root.

While the function definition may appear complicated, much of the complication is to ensure that we do not go off into some loop.  The basic idea is quite simple:  

First find the approximate root using FindRoot, calculating the step size between iterations as it works.  

Then, add the number of extra approximate digits of accuracy to the AccuracyGoal option and call FindRoot again, repating this last step until the step size criterion is satisfied.  

The number of extra digits need for the AccuracyGoal is approximated by first approximating the order of the root (p = Log[change in function]/Log[change in x]) and multiplying this by the difference of the original AccuracyGoal and the Log base 10 of the current step.

StepErrorFindRoot definition

StepErrorFindRoot[f_,{x_,x0_},opts___] :=
Module[{},
	wp = WorkingPrecision /. {opts} /. Options[FindRoot];
	wp = Max[wp,$MachinePrecision];
	ag = AccuracyGoal /. {opts} /. Options[FindRoot];
	If[ag === Automatic,ag = wp-10];
	xold = x0;
	step =.;
	Print["Using FindRoot with AccuracyGoal->",ag,
		", WorkingPrecision->",wp];
	sol = FindRoot[fnew = f;
		If[NumberQ[x] && (x != xold),
			step = Abs[x-xold];
			fstep = Abs[fnew-fold];
			xold = x];
		fold = fnew,{x,x0},opts];	
	trys = 0;
	agf = ag;
	While[((trys++) < 20) && NumberQ[step*fstep] && step > 10^(-ag),
		p = Max[1,Log[fstep]/Log[step]];
		agf = agf + Ceiling[p*(Log[10,step]+ag)];
		wp = Max[$MachinePrecision,wp,agf*2];
		xnew = x /. sol;
		step =.;
		Print["Using FindRoot with AccuracyGoal->",agf,
		", WorkingPrecision->",wp];
		xx = x /. sol;
		sol = FindRoot[fnew = f;
			If[NumberQ[x] && (x != xold),
				step = Abs[x-xold];
				fstep = Abs[fnew-fold];
				xold = x];
			fold = fnew,{x,xx},
			AccuracyGoal->agf,
			WorkingPrecision->wp,
			DampingFactor->p,opts]];
	If[!NumberQ[step*fstep],
		If[Abs[x - xx] > 10^(-ag),
			Print["Non numeric step: 
			failed to meet the step size accuracy goal"]],
	 	If[step > 10^(-ag),
	 		Print["Failed to meet step size accuracy goal"]]];
	sol]

mypi  = x /. StepErrorFindRoot[Cos[x]+1,{x,3}]

N[Pi - mypi]

Which is certainly within the expected tolerance of 10^-6.

PrecisionGoal vs AccuracyGoal: NIntegrate

NIntegrate works by dividing the interval of integration into subintervals.  On each subinterval an integration approximation is used.  Associated with the integration approximation on each interval is an error estimate. The interval error estimates are added up to get the global error estimate.   If the global error exceeds the error tolerance, intervals are subdivided, starting with the ones with the largest error until the tolerance is met.

The following set of examples highlights the difference between specifying an AccuracyGoal and specifying a PrecisionGoal.  The difference is clear when the size of the integrand is rescaled.

err[b_,opts___] := (c = 0;
		NIntegrate[c++;b Sqrt[x],{x,0,100},opts] - 
		Integrate[b Sqrt[x],{x,0,100}]);
	B = Table[10^j,{j,0,10}];

The default is an AccuracyGoal of Infinity, so that the error control is based on the PrecisionGoal which is by default 6 for machine precision calculations.  (The variable c is a counter which keeps track of the number of times NIntegrate needs to evaluate the function)

TableForm[{#,err[#],c}& /@ B,
	TableHeadings->{None,{"b","error","c"}}]

It is clear that as the size of the integrand increases, the error increases correspondingly.  In fact, the function is evaluated at the same points and the only difference between the results is the magnitude of the solution.

On the other hand, when the error estimate is based on AccuracyGoal, the method tries to achieve the same absolute error tolerance independent of the size of the intgrand.

TableForm[
	{#,err[#,PrecisionGoal->Infinity,AccuracyGoal->6],c}& /@ B,
	TableHeadings->{None,{"b","error","c"}}]

The method suceeds in meeting its accuracy goal until it gets to the point where machine precision numbers are not sufficient to do the calculation to the desired accuracy. Redoing the calculation with a higher WorkingPrecision will eliminate this problem.

Local vs Global Error Estimates: NIntegrate and NDSolve

Consider a simple integration problem.  Let

Clear[f];
f[t_] := Sin[t]

We will focus on two different numerical methods to obtain an approximation to the value of the integral,

xe[t_] = 2+Integrate[f[t],t]

as a function of the variable t.  Note that the integration constant was chosen so that NDSolve could be given the same values for AccuracyGoal and PrecisionGoal as NIntegrate has for its defaults (Infinity and 6).  NDSolve cannot meet the default PrecisionGoal of 6 if the initial condition is 0.

Naturally, the first method that comes to mind for approximating the value of an integral is NIntegrate.  A function which finds an approximation as a function of t can easily be defined as follows:

Off[NIntegrate::nlim];
xni[t_] = NIntegrate[f[tv],{tv,0,t}]+1

However, this function has the significant drawback that for each value of t the function is evaluated at, an approximation of the integral with NIntegrate is necessary.  If more than a few evaluations are necessary (say to make a plot), the time becomes a substantial problem.

On the other hand, this is also a special case of the class of problems that NDSolve is designed to solve, and NDSolve returns an InterpolatingFunction object which approximates the value for a given t by the process of interpolating the data collected from NDSolve, so this is evaluated quickly for each value of t.

sol = NDSolve[{x'[t] == f[t],x[0] == 1.},
	x,{t,0,160*Pi},
	PrecisionGoal->6,
	AccuracyGoal->Infinity,
	MaxSteps->5000];
xn = x /. sol[[1]]

The versatility gained by using NDSolve--an approximate result with values at any t versus the need to compute for each value of t needed--does not come without price.  For smooth problems, the major loss is that of global error control: while NIntegrate surveys the error over the entire interval, NDSolve only considers the error in one part of the solution at a time.  In other words, NDSolve implements only local error control.

The following table demonstrates the difference between the two algorithms:

Off[NIntegrate::ncvb];
TableForm[tab = Table[
		{t,N[xni[t]-xe[t]],N[xn[t]-xe[t]]},
		{t,20*Pi,160*Pi,20*Pi}],
	TableHeadings->{None,{"t","NIntegrate\nerror","NDSolve\nerror"},None}]

It is evident from the table that the error from the NDSolve approximation tends to grow as the interval length (and hence the number of steps) grows.  Each step is taken so that a (generally conservative) estimate of the error doesn't exceed the default PrecisionGoal of 6.  However, with many steps, it is quite possible that the total error does exceed the default PrecisionGoal.  The NIntegrate approximation has error which varies, but doesn't show such a tendancy to grow as the interval is increased in length.  There is of course a limit to how long the interval can get before the NIntegrate error begins to grow--there are factors such as round-off error and limitations on the number of bisections done by NIntegrate to achieve a result which will become important for extremely long intervals.

eps = 10^(-$MachinePrecision);
ttab = Transpose[tab];
Show[
	LinearLogListPlot[
		Transpose[{ttab[[1]],eps+Abs[ttab[[2]]]}],
		PlotJoined->True,
		DisplayFunction->Identity],
	LinearLogListPlot[
		Transpose[{ttab[[1]],eps+Abs[ttab[[3]]]}],
		PlotJoined->True,
		PlotStyle->Red,
		DisplayFunction->Identity],
	PlotRange->All,DisplayFunction->$DisplayFunction];

The error from NDSolve is will not always grow like this.  In fact, it is more typical that errors will cancel out and the growth will be mcuh smaller than this.  It is important to realize, however, that the error can, and sometimes does, grow in this fashion.

Sensitivity to Initial Conditions: NDSolve

Some differential equations can be very sensitive to initial conditions.  For equations where this is possible, you may want to check your solution carefully.

For example, find a solution to the pendulum equation.

sol = NDSolve[{x''[t] +Sin[x[t]] == 0,
	x[0] == 0, x'[0] == 2},x,{t,0,50}][[1,1,2]]

This plot seems a bit strange:  you expect the undamped and  unforced pendulum equation to produce solutions which either oscillate or grow (the case where the pendulum has enough energy to swing in a circle--the angle grows), but not to do some of both as appears the case in the plot below.

Plot[sol[t],{t,0,50}];

Since this solution is a bit suspect, it is worth checking to see what happens when the initial condition is perturbed slightly.  Here the size of the perturbation is similar to the size of the local truncation error NDSolve uses based on the default values of AccuracyGoal and PrecisionGoal.

psol = NDSolve[{x''[t] +Sin[x[t]] == 0,
	x[0] == 0, x'[0] == 2-10^-6},x,{t,0,50}][[1,1,2]]

Plot[{sol[t],psol[t]},{t,0,50},PlotStyle->{Red,Blue}];

It is clear that the two solutions are quite different even though they started very close together.

What has happened is that the initial condition lies on the separatrix for the pendulum equation which separates the region of phase space where oscillation occurs and the region where the pendulum swings in circles.  Along this separatrix is the (unstable) equilibrium point where the pendulum is vertically upwards--any perturbation from this point and the pendulum swings.  The initial condition given is equivalent to giveng just enough velocity so that the pendulum will rise to the top and asymptotically approach the equilibrium point.  Physically, any small perturbation would eventually move the pendulum away from the equilibrium point.  In solving the equation approximately, we have introduced perturbations which are of size roughly 10^-6--the size of the error allowed by the default accuracy and precision goals.

One way to reduce these perturbations is to increase the default values of AccuracyGoal and PrecisionGoal.  It is not possible to increase them very much without using higher precision numbers. This can be done by setting the WorkingPrecision option.  When the WorkingPrecision is greater than machine precision, AccuracyGoal and PrecisionGoal Automatic give default values of 10 minus the WorkingPrecision.  Thus, the command

asol = NDSolve[{x''[t] +Sin[x[t]] == 0,
	x[0] == 0, x'[0] == 2},x,{t,0,50},
	MaxSteps->Infinity,
	AccuracyGoal->22,
	PrecisionGoal->22,
	WorkingPrecision->32][[1,1,2]];

computes a solution which should have the local truncation error near 10^-22.

As can be seen from the plot, the solution stays near the real one for much longer before it deviates.

Plot[asol[t],{t,0,50},PlotRange->All];

In reality, there is no way to compute an approximate solution with NDSolve which will be accurate for all time.

Exact Solution

exsol[t_] = x /. Solve[t == EllipticF[x/2,1],x][[1]]

There are many situations where error can be greatly magnified due to sensitivity to initial conditions.  If you think that a problem you are solving may have this property, it is a good idea to see what perturbed solutions do, or do other analysis to check on the accuracy of the solution.  In some cases, going to higher precision will help, but not always.