Numerical Mathematica

Main Theme:
Using Accuracy and Precision

Mathematica's numerical commands allow excellent control over result quality and usually the Automatic settings are sufficient for most uses.  One situation which arises where the Automatic settings are not ideal is when you are using more than one numerical command at a time -- e.g. using one to provide a function value for another.   For these types of applications to be efficient, you need to have a good understanding of the meaning of now Mathematica interprests accuracy and precision and how these are related to the AccuracyGoal and PrecisionGoal options.  In this talk, I will guide you first through an introduction to accuracy and precision in Mathematica and then two examples where the concepts are used for tight error control.

Initialization

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

$MinPrecision  =  -Infinity;

$Path  =  
        Join[$Path,  {"C:\\Documents and Settings\\rknapp\\My \
Documents\\Talks\\2001 IMS", "C:\\Documents and Settings\\rknapp\\My \
Documents\\Talks\\2001 IMS\\Packages"}];

Needs["Utilities`Scale`"]

Basics of Accuracy and Precision

A starting point is to look at the relevant section from The Mathematica Book: Section 3.1.4.

The concept becomes clearer if you use the concept of scale to form a relationship between the two.
The scale of a number is the number of digits it has to the left of the decimal point.  Thus, for any number on the real line x, we have

Precision[x] == Scale[x] + Accuracy[x]

Here is a pictorial representation of the quantities:

(Overscript[(Overscript[(Underscript[\(\(\(\ \)\(x\_1\ x\_2\ ⋯\ x\_s\)\)\+⏟\), Scale]).(Underscript[\(\(x\_\(s + 1\)\ x\_\(s + 2\)\ ⋯\ x\_\(s + a\)\)\+⏟\), Accuracy]), ⏞]), Precision])

Any of the values for Accuracy, Precision and Scale can be negative.

Lets see how this relationship holds for some examples

PSA[x_]  :=  {Precision[x],  Scale[x],  Accuracy[x],  
        Precision[x]  == (Scale[x]  +   Accuracy[x])}

The Mathematica Book example ...

PSA[xacc = 431.123145333555141444]

... almost works.  The problem is that you really need to consider Precision and Accuracy as approximate quantities, not integral numbers of digits:

PSA[x_]  :=  {Precision[x,  Round -> False],  Scale[x],  
        Accuracy[x,  Round -> False],  
        Precision[x,  
            Round -> 
              False]  == (Scale[x]  +   Accuracy[x,  
                Round -> False])}

and now...

PSA[xacc]

... comes out as it should!

This even works for exact numbers

PSA[1]

Notice that for small (<< 1) numbers, the scale is negative

PSA[10. ^-1000]

Since it clearly makes more sense to think of Precision and Accuracy as floating point quantities than it does as a particular number of digits.  A  number with a particular accuracy can best be thought of as representing an interval which is centered at the number.

Suppose that x has Precision[x] == p  and Accuracy[x] == a.  Then the absolute width of the error interval for x is (10^(-a)).  This means that x represents any number \(TraditionalForm\`\(x\_*\)\)in the interval

\(TraditionalForm\`x - \(1\/2\) 
          10\^\(-a\) ≤ \(x\_*\) ≤ 
      x + \(1\/2\) 10\^\(-a\)\)

or equivalently

\(TraditionalForm\`x(
        1 - \(1\/2\) 10\^\(-p\)) ≤ \(x\_*\) ≤ 
      x(1 + \(1\/2\) 10\^\(-p\))\)

To make this concrete, consider the number xacc defined above, where the error interval has width

10^-Accuracy[xacc,  Round -> False]

and so xacc could be any number \(TraditionalForm\`\(x\_*\)\) in

\(\*\(xinf\  = \ 
        SetPrecision[xacc, \ Infinity];\)\)
\(\*\(iinf\  = \ 
        SetPrecision[10^\(-Accuracy[xacc]\), \ Infinity]/
          2;\)\)
\(\*xinf\  - iinf\  ≤ \ \(x\_*\) ≤ \ 
      xinf\  + \ iinf\)

Rational numbers have been used here to make the boundary of the interval sharp and clear.

Here is an example which illuminates the concept of negative precision.  Consider

poorx  =  3`-5

This could be any number in the interval

\(\*\(pinf\  = \ 
        SetPrecision[SetPrecision[poorx, \ 100], \ 
          Infinity];\)\)
\(\*\(pint\  = 10^\(-Precision[poorx]\);\)\)
\(\*pinf*\((1\  - pint/2)\)\  ≤ \ \(x\_*\) ≤ \ 
      pinf*\((1\  + \ pint/2)\)\)

The error interval size is consistent with the accuracy

150003  +  149997  ==  
      10^-Accuracy[poorx,  Round -> False]

While the idea of precision breaks down around zero, accuracy is absolute so it can be used anywhere.

If you would like to learn more about these ideas and how they are used to form our significance arithmetic, and excellent resource is Mark Sofroniou's talk on Efficient and Numerically Stable Computations from the 1998 Chicago conference.

PrecisionGoal and AccuracyGoal

The concepts of accuracy and precision are extended to numerical functions with options AccuracyGoal and PrecisionGoal.  An example helps illuminate this

NIntegrate[Sin[2 π x], {x, 0, 1}]

Gives a warning message because it could not meet its default criteria:

Options[NIntegrate]

The result may seem quite good, but in terms of relative error it is poor.  The appropriate criterion to use here is absolute error, so you could use

NIntegrate[Sin[2 π x], {x, 0, 1},  
      AccuracyGoal -> 6]

And we are fine.  In general, when both Precision and AccuracyGoal are present, the result is deemed OK if the least stringent of the two is estimated to be met.

In general, Precision is the better criterion, since it allows both

NIntegrate[Exp[x], {x, 100,  101} ]

and

NIntegrate[Exp[x], {x, -101,  -100}]

To give reasonable answers with the same relative (though certainly not absolute!) error.

Simple Example: Division

Simple division can be done using multiplication and addition only by using Newton's method.  With

f[a_][x_]  :=  a  -  (1/x)

then

Expand[x  -  f[a][x]/f[a]'[x]]

so you can define

DivisionStep[a_][x_]  := (2  - a x) x

Then you can look for the fixed point

pin  =  N[Pi];
FixedPoint[DivisionStep[pin],   .1]

Machine numbers were used in the computation above because it was started with a machine number.  If you start with a bignum, then you can get a higher precision result.

prec  =  100000;
pin  =  N[Pi,  prec];
start  =  N[1/10, prec];
Timing[rip  =  
          FixedPoint[DivisionStep[pin],  start];]

which takes a little while, but does give a result which is close

rip*pin  -  1

Newton's method is a self correcting algorithm in the sense that if the iterations are converging, errors will disappear.  Thus, it is possible to start with lower precision and increase the precision as you get closer to the root.

AdaptiveNewton[step_,  start_,  
        opts___]  :=  
Module[{goal,  
          prec}, 
c  =  0; 
        goal  =  
          PrecisionGoal  /.  
            Flatten[{opts,  
                Options[AdaptiveNewton]}]; 
prec  =  
          Min[Precision[start],  goal]; 
SetPrecision[
          NestWhile[AdaptiveNewtonStep[step,  prec],  start,  
            Function[prec  <  2*goal]],  
          goal]
]

FixedPoint could be used, but it is not ideal, because its criterion for stopping is that the result is no longer changing up to the least common precision.  Since the precision is being changed, it would not terminate.  NestWhile fits the bill very nicely because it allows you to put in a function to test for a stopping condition.  Using it, we need to define a function which handles each iterative step.

\(\(AdaptiveNewtonStep[step_, \ prec_]\)[\ x_]\  := \ 
      Module[{xnew, \ Δx}, 
xnew\  = \ 
          step[x]; 
Δx\  = \ 
          xnew\  - \ 
            x; 
If[\((Δx\  == \ 
                0)\) || \ \((prec\  < \ 
                2 \((Scale[
                        x] - \ \ Scale[Δx])\))\), 
\ 
          prec\  *= \ 
            2. ; 
print["Increasing Precision to ", \ 
            prec, \ " at iteration ", \ \(++c\)]; 
xnew\  = \ SetPrecision[xnew, \ prec]]; 
xnew]\)

There two useful tricks being used here.

Arguments by reference

First, we are changing prec with TimesBy.  This requires that prec be a symbol (and also that it has a value, since if not a recursion error will occur).  But if you try it

prec  =  $MachinePrecision; {AdaptiveNewtonStep[
          DivisionStep[pin],  prec][1. /Pi], prec}

That is because by default arguments are evaluated.  You might expect that it would be possible to use Unevaluated, but this gives

{AdaptiveNewtonStep[DivisionStep[pin],  Unevaluated[prec]][
        1. /Pi], prec}

Using Unevaluated can be difficult at best.  Since any use of the function requires that prec not be evaluated, it is better to associate this with the function itself anyway.

SetAttributes[AdaptiveNewtonStep,  HoldAll];

Now everything works just as we want.

{AdaptiveNewtonStep[DivisionStep[pin],  prec][1. /Pi], 
      prec}

Even if the step were a complicated function which should be evaluated first, that is not a problem since it will have already been evaluated when using AdaptiveNewton

In Mathematica, the default is to pass arguments by value.  For example

Subtract[Times[rip,  Pi],  1]

is computed by descending the expression recursively, evaluating arguments, then evaluating the functions.  The code for Times and Subtract have no (nor do they need to) knowledge about whether rip, Pi, and 1 are associated with variables.

On the other hand,

Set[a,  5]

requires that a be a symbol.  a could have a previous value, which will be replaced by 5.  Set works with the symbol a, which after set is run is a reference to the number 5.  In Mathematica references are associated with symbols.

What about

b[1]  =  a

In this case it appears that we are assigning to b[1], which is not a symbol.  In reality, the rule is attached to b:

DownValues[b]

However, it is important to note that it is possible to pass by reference in Mathematica by delaying evaluation and this is a useful technique to be aware of.

Using approximate scale

The inequality,

prec  <  2 (Scale[x] -   Scale[Δx])

probably seems quite obscure, but if you use the definition of scale,

prec  <  
      2 (Log[10. ,  Abs[N[x]]]  -   Log[10. ,  
              Abs[N[Δx]]])

and you divide both sides of the equation by -2 and exponentiate, this becomes

Abs[Δx]/( Abs[x]) <  
      (10^(-prec/2))

The left-hand side is the size of Δx relative to x and the right hand side represents an (accuracy) tolerance of half of the digits available in prec.

Why half?  This is tied in with our doubling of the precision, and is associated directly with the quadratic convergence of Newton's method: If (Subscript[x, 1]), (Subscript[x, 2]), ...  are a sequence of Newton iterates converging to a root limit \(TraditionalForm\`\(x\_*\)\), then

\[LeftBracketingBar](Subscript[x, k + 1]) - 
          (Subscript[x, k])\[RightBracketingBar] -> \[LeftBracketingBar](Subscript[x, k]) - (Subscript[x, k - 1])(\[RightBracketingBar]^2)

Thus, if Δ(Subscript[x, k]) = |(Subscript[x, k]) - 
          (Subscript[x, k - 1])| is of order (10^(-prec/2))(i.e. the first prec/2 digits cancel in the difference) you expect the next Δ(Subscript[x, k + 1])= |(Subscript[x, k + 1]) - 
          (Subscript[x, k])| to be of the order of (10^(-prec)), so if you want (Subscript[x, k + 1])to have any significant digits different than (Subscript[x, k]), you had better increase the precision!

The idea of using the scale is that it avoids doing any extra computations at the full current precision.  For example, doing the division in the derived inequality above would be as expensive as the problem we're trying to solve!

Note that it is important to keep in mind, that using scale is only appropriate where an approximate cutoff is acceptable.  If you need a condition to satisfied strictly, you should check that condition at full precision.

Now lets try it to see if it works

Timing[
      rip  =  AdaptiveNewton[DivisionStep[pin],  N[1/Pi],  
          PrecisionGoal -> 100000]; c]

It turns out we can do better yet.  Notice that the final precision used exceeds the goal precision.  This means we have used too much precision.  Since multiplication cost increases quickly with precision, this will affect the overall timing.  A solution is to choose the sequence of precision carefully so we wind up with exactly the goal precision.

AdaptiveNewton[step_,  start_,  
        opts___]  :=  
Module[{goal,  
          prec}, 
c  =  0; 
        goal  =  
          PrecisionGoal  /.  
            Flatten[{opts,  
                Options[AdaptiveNewton]}]; 
prec  =  
          Min[Precision[start],  16]; 
prec  =  
          2^(Floor[Log[2. ,  prec]]  +  
                FractionalPart[
                  Log[2. ,  goal]]); 
SetPrecision[
          NestWhile[AdaptiveNewtonStep[step,  prec],  start,  
            Function[prec  <  2*goal]],  
          goal]
]

Doing this saves one extra evaluation, and since it was the most expensive one, a substantial proportion of the time is saved.

Timing[
      rip  =  AdaptiveNewton[DivisionStep[pin],  N[1/Pi],  
          PrecisionGoal -> 100000]; c]

This compares favorably with Mathematica's internal division algorithm

Timing[1/pin;]

It is worth noting that this is only a small multiple of the cost of a multiplication at this precision:

Timing[rip*rip;]

While this is quite a simple calculation, you have seen that some extra thought about controlling precision can bring dramatically better performance.   In more complicated calculations, the control of precision and accuracy is often much less clear, but can be just as important.

Another important conclusion is that it is also possible to write general code without sacrificing speed.  The code above is almost as fast as Mathematica's internal division code, but it is not special to division.  For example, we could use it to find SquareRoots.

SquareRootStep[a_][
        x_]  :=  (x^2  +  a)/(2 x)

Timing[rip  =  
          AdaptiveNewton[SquareRootStep[2], Sqrt[2. ],  
            PrecisionGoal -> 100000];]

but it is quite slow compared to the internal algorithm

Timing[N[Sqrt[2],  100000];]

This, was expensive because each iteration included a division.  Using the function

a - (1/(x^2))

for the reciprocal of the square root, it is possible to construct a divisionless step (division by 2 is simple for Mathematica because with our binary representation, it is just a bit shift.

ReciprocalSquareRootStep[a_][x_]  :=  
      x*(3  -  a x^2)/2

Using the divisionless step, it is not a lot slower than the internal algorithm (the internal one avoids the division at the end by taking a special last step)

Timing[rip  =  
          AdaptiveNewton[ReciprocalSquareRootStep[2], Sqrt[2. ],  
            PrecisionGoal -> 100000];]

More involved example: Cumulative integrals.

This is a question which comes up from time to time.  If you have a given function, f, how do you compute

(∫_a^x) f(
            s) \[DifferentialD]s ==  c

numerically a value of x so that the equation is satisfied.

Lets take a simple example

f[x_]  :=  (1/(Sqrt[π])) Exp[-x^2]

with cumulative integral

Cf[x_,  opts___]  :=  
      NIntegrate[f[s], {s, -Infinity,  x},  opts]

and look for what x limit gives 99.9% of the full value of 1 over the real line:

q  =  999/1000;

graphically, this looks like

Plot[
        Evaluate[{Integrate[f[s], {s, -Infinity, x}], q}], {x, -3, 
          3}];

The simplest way to do it is with FindRoot:

root  =  
      x  /. FindRoot[Cf[x]  -  q, {x, 0},  
          Jacobian -> f[x]]

Cf[root]  -  q

Now suppose we want to know x to p digits of precision.

p  =  10;

Presumably, all we have to do is change the AccuracyGoal option in FindRoot to AccuracyGoal->p

root  =  
      x  /.  FindRoot[
          NIntegrate[f[s], {s,  -Infinity,  x}]  -  q, {x, 0},  
          Jacobian -> f[x],  AccuracyGoal -> p  +  2]

However, a plot shows that our solution does not lie in the requested error interval.

Plot[
        Cf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p}];

Experimentally, you can find that increasing the AccuracyGoal to AccuracyGoal->p + 2 takes care of this problem:

root  =  
        x  /.  FindRoot[
            NIntegrate[f[s], {s,  -Infinity,  x}]  -  q, {x, 0},  
            Jacobian -> f[x],  
            AccuracyGoal -> p  +  2];
Plot[
        Cf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p}];

Now compare it to the exact solution. (Yes, this example was conveniently chosen to have an exact solution.  Such a solution is often very useful when trying to work out the details of an algorithm.)

ExactCf[x_]  :=  (Erf[x]  +  1)/2

Plot[
        ExactCf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p},  PlotRange -> All];

Ouch!  The plot shows that we are NOT within the requested error interval for the exact solution.  In fact, an expanded plot shows that our error is a factor of 10 too large.

Plot[
        ExactCf[(1  +  x) root] - q, {x, -10^(1 - p), 
          10^(1 - p)}];

The entire idea of a numerical command like this is to use error estimates conservative enough to be fairly sure we are within the prescribed error interval from the exact solution.

Clearly something is wrong with the approach taken so far.  What is it?

This plot demonstrates the source of the problem quite clearly:

Plot[
        Cf[(1  +  x) root]  -  
          ExactCf[(1  +  x) root], {x, -10^(1  -  p),  
          10^(1  -  p)},  PlotRange -> All];

The value of the function FIndRoot is trying to find a root for is varying in essentially random way on the same scale as the size of the function.  There is no way FindRoot can be expected to find the correct value under these circumstances.

If you ask NIntegrate to do a better job, then we can reduce the scale of local variation:

Plot[
        Cf[(1  +  x) root,  AccuracyGoal -> 12,  
            PrecisionGoal -> 12]  -  
          ExactCf[(1  +  x) root], {x, -10^(1  -  p),  
          10^(1  -  p)},  PlotRange -> All];

This should be sufficient to get the root to the desired precision

root  =  
      x  /.  FindRoot[
          Cf[x,  AccuracyGoal -> 12,  PrecisionGoal -> 12]  -  
            q, {x, 0},  Jacobian -> f[x],  
          AccuracyGoal -> p +  2]

and a plot shows that it was.

Plot[
        ExactCf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p},  PlotRange -> All];

So, for this problem, we have a balancing act between the AccuracyGoal and PrecisionGoal for NIntegrate, the AccuracyGoal for FindRoot, and the precision or accuracy we need for the root.  Because the cost of NIntegrate goes up exponentially fast with increase Accuracy/Precision Goals, it is important to get parameters quite close to the balance point.

When you have an exact result, it is not terribly difficult to use graphs to determine the constraints needed to obtain an adequate solution.  It is a more challenging problem to determine the constraints automatically with no exact solution to compare with.  Below, I develop a general (though not completely bulletproof) way of handling the balance of parameters.

The local sensitivity around the root is given by the slope.  The local linear model around the root \(TraditionalForm\`\(x\_*\)\) is f(x) ≅ (x - \(TraditionalForm\`\(x\_*\)\)) \(TraditionalForm\`f' \((\(x\_*\))\)\), and so the error in x is approximately x - \(TraditionalForm\`\(x\_*\)\) ≅ \(TraditionalForm\`\(f(x)\)/
            f' \((\(x\_*\))\)\  ≅ \ \(f(x)\)/f' \((x)\)\).  The last ratio should, of course, look very familiar since it is plain and simple the Newton step!

If you write your own Newton iteration loop, you could use an appropriate bound on this as a stopping condition (much as we did for the division example).  Currently, since FindRoot only uses the size of the function as a stopping condition, you need to adjust the AccuracyGoal according to the slope.  Typically a reasonable approximation to the root will give you a good enough approximation to the slope at the root to adjust the option.  This can be done with a two step operation.

crude  =  
        x  /.  FindRoot[Cf[x]  -  q,  {x, 0},  
            Jacobian -> f[x]];
slope  =  f[crude]

Given that, we want the ratio of the function value to the slope to be less than the precision we are aiming for.  In particular

((10^(-FindRootAG))/slope) ≤ 
      x (10^(-PG))

The minimal condition is at equality and taking logarithms on both sides, we have

-FindRootAG - (Subscript[log, 10])(slope) == (Subscript[log, 10])(
          x) - PG

Recalling the definition of Scale, you get

FindRootAG  =  p  -  Scale[slope]  -  Scale[crude]

Before calling FindRoot again, it is appropriate to think about what we have to require from NIntegrate.  This is actually relatively simple:  we are asking FindRoot to iterate until the difference of the value of the integral and  a constant be of scale equal to -FindRootAG.  Assuming the constant has sufficient precision (if it doesn't you're sunk in any case!), then that is the minimal accuracy which is required from NIntegrate.  Actually, to be sure that the function appears smooth (and not random like the plots above) to FindRoot near the approximate root, it is a good idea to give a little extra.  This is also helpful to put us on the safe side since, of course, NIntegrate is basing its error control on estimates as well.

NIntegrateExtraAccuracy  =  1;

With all this, we can compute our final root with

root  =  
      x  /.  FindRoot[
          Cf[x,  AccuracyGoal -> 
                FindRootAG  +  NIntegrateExtraAccuracy,  
              PrecisionGoal -> Infinity]  -  q,  {x,  crude},  
          Jacobian -> f[x],  
          AccuracyGoal -> FindRootAG]

And from the plot, it looks like we got it!

Plot[
        ExactCf[(1  +  x) root]  -  q, {x, -10^-p, 
          10^-p},  PlotRange -> All];

Now that we have done all of this head scratching, it's worth putting together a function which combines these ideas.  I have put together a package which merges the commands,

FindIntegralLimit1

Needs["FindIntegralLimit1`"]

?FindIntegralLimit1

A utility this package uses which is worth discussing is encapsulated in the package AutomaticPGAG.m, which mimics the defaults for options AccuracyGoal->Automatic and PrecisionGoal->Automatic which Mathematica's built-in numerical functions use.

An advantage of a nice simple function is that it is easy to use it in operations, such as plotting:

Plot[
        x  /.  FindIntegralLimit1[f[s], {s, -Infinity,  x,  0},  
            q],  {q, 1/2, 9999/10000}];

Shows the expected behavior.  What is more interesting is a plot of the error:

Plot[
        Log[10, Abs[
            InverseErf[0,  
                2 q  - 1]  -  (x  /.  
                  FindIntegralLimit1[f[s], {s, -Infinity,  x,  0},  
                    q])]],  {q, 1/2, 9999/10000},  
        PlotRange -> All];

The structure of the plot is an artifact of the way error control is handled.  It is not directly a feature of the problem being solved.

You may have noticed that this took a while to produce.  The speed of the command becomes more of an issue when more precision is requested.

Timing[
      FindIntegralLimit1[f[s], {s, -Infinity,  x,  0},  999/1000,  
        WorkingPrecision -> 24]]

One way this can be made faster is to rewrite the input so that NIntegrate has less work to do:

Timing[
      FindIntegralLimit1[f[s], {s, Infinity,  x, 2.185124219133004},  
        999/1000  -  1,  WorkingPrecision -> 24]]

but that is not always possible.  The next section addresses some possibilities for doing this.

Some possible speedups

Cache computed integral values

It is clear that doing the smaller integral was much faster for NIntegrate.  By limiting the extent of the region over which NIntegrate has to work, we should be able to save some time.  Because

(∫_a^x) f(
            s) \[DifferentialD]s == (∫_a^((Subscript[x, 1])))f(
              s) \[DifferentialD]s + (∫_((Subscript[x, 1]))^x) f(
              s) \[DifferentialD]s

We can do the following

cachex  =  -Infinity;
cachei  =  0;
Cfcache[x_,  
          opts___]  :=  (cachei  +=  
            NIntegrate[f[s], {s, cachex,  x},  opts];  cachex  =  x;  
          cachei);

Timing[Cf[0,  WorkingPrecision -> 24]]

Timing[Cfcache[0,  WorkingPrecision -> 24]]

Timing[
      Cf[2.185124219133004264123974081742574`24,  
        WorkingPrecision -> 24]]

Timing[
      Cfcache[2.185124219133004264123974081742574`24,  
        WorkingPrecision -> 24]]

You have to be careful with this approach, however because it is not quite as accurate.  NIntegrate computes a global error estimate.  This is lost  when you add the results of several invocations in NIntegrate.   However, if the number added together is not too large, the extra error is not too large.  For example, even if you add 10 intervals together, we can expect the same quality result by adding 1 extra to NIntegrateExtraPrecision.

I've put this into a second version of the package, FindIntegralLimit2.m.

Needs["FindIntegralLimit2`"]

and, sure enough, it is faster on this example.

Timing[
      FindIntegralLimit2[f[s], {s, -Infinity,  x,  0},  999/1000,  
        WorkingPrecision -> 24]]

Use NDSolve

NDSolve can be used as a numerical integrator also, though since it does not have global error control, it is usually not appropriate.  For the crude approximation, however, it can in some cases work quite nicely.

Since NDSolve cannot handle infinite limits, there are two things we can do.  One is to do a variable transformation much like what NIntegrate does automatically.  While this is not exceptionally difficult, I will not cover this method here.

The easiest is to identify a finite part of the integral to use.  So that you can see how NDSolve works, lets change the problem so the lower limit is 0.

ndsol  =  
        First[u  /.  
            NDSolve[{u'[s]  ==  f[s],  u[0]  ==  0},  
              u,  {s, 0,  10}]];

Plot[ndsol[s], {s, 0, 10},  PlotRange -> All];

FindRoot[ndsol[s]  -  499/1000, {s, 1, 0, 10}]

To make this a little more efficient, NDSolve has a neat little option called StoppingTest, which causes NDSolve to stop when a given condition evaluates to True.

ndsol  =  
      First[u  /.  
          NDSolve[{u'[s]  ==  f[s],  u[0]  ==  0},  
            u,  {s, 0,  10},  
            StoppingTest -> (u[s]  >  499/1000)]]

This means that in the last step to s == 2.32251, u crossed the threshold.  So ideally, we want to narrow the FindRoot search to between

ndend  =  ndsol[[3, 1, -1]]

and

ndprev  =  ndsol[[3, 1, -2]]

When FindRoot is called, it will use a variant of Brent's algorithm which keeps the root bracketed.

FindRoot[ndsol[s]  -  499/1000,  {s,  ndprev,  ndend}]

Because of the fact that NDSolve does not have global error control, we do not want to use it for the final approximation.  However, this is quite adequate for the crude approximation.  The use of NDSolve has been incorporated into a third version of the package, FindIntegralLimit3.m

Needs["FindIntegralLimit3`"]

Timing[
      FindIntegralLimit3[f[s], {s, -Infinity,  x,  0, 10},  
        999/1000,  WorkingPrecision -> 24]]

This approach could save more time for an example where the initial root isolation was more difficult.

FindIntegralLimit

So what should the final package have that I have not included in any of my three versions so far.

• It should have a tidy way of switching methods for finding the crude approximation.  One possibility is suggested by the extra argument required by FindIntegralLimit3.  A drawback of this is that simply adding an extra argument is not intuitively an obvious way to change to using NDSolve.

• When the precision of the starting condition is bigger than machine precision, the crude approximation should likely be done with that precision, similar to the way FindRoot works.  That way, if you have a function which gives garbage for machine numbers, you are not wasting time with useless calculations.

• It should check carefully for failures in NIntegrate, NDSolve  and FindRoot and exit gracefully with a tasteful message.

• When the precision will have to be increased several times in the last computation, it would be good to have an automatic way to have to cache integral restarted whenever the precision is changed.    (In versions 2, 3 we go to full precision for the last call to FindRoot to avoid the problem of FindRoot increasing precision.  Such a large increase in precision will be rare since NIntegrate limits the practical precision which can be used.

• It would be good to come up with alternative schemes in the case where the slope (Jacobian) is zero or nearly so at the crude approximation.

• Tests.  I have just shown the idea of the function on a relatively simple example.  To get a robust package working, you would have to come up with a variety of tests.

Summary

Once you have a good understanding of the concepts of accuracy and precision, you can gain a good understanding of how good you should expect solutions from numerical functions to be when you specify values for the AccuracyGoal or PrecisionGoal options.  From this, it is possible to write more complex numerical functions which combine several mroe elementary numerical operations.