The `Richardson`

package provides a function `extrapolate`

that
extrapolates any function `f(x)`

to `f(x0)`

, evaluating
`f`

only at a geometric sequence of points `> x0`

(or optionally `< x0`

) or at a given sequence of points. `f(x)`

can
return scalars, vectors, or any type implementing a normed vector space.

The key algorithm is Richardson extrapolation using a Neville–Aitken
tableau, which adaptively increases the degree of an extrapolation
polynomial until convergence is achieved to a desired tolerance
(or convergence stalls due to e.g. floating-point errors). This
allows one to obtain `f(x0)`

to high-order accuracy, assuming
that `f(x0+h)`

has a Taylor series or some other power
series in `h`

. (See e.g. these course notes by Prof. Flaherty at RPI.)

```
extrapolate(f, h; contract=0.125, x0=zero(h), power=1,
atol=0, rtol=atol>0 ? 0 : sqrt(ε), maxeval=typemax(Int), breaktol=2)
```

Extrapolate `f(x)`

to `f₀ ≈ f(x0)`

, evaluating `f`

only at `x > x0`

points
(or `x < x0`

if `h < 0`

) using Richardson extrapolation starting at
`x=x₀+h`

. It returns a tuple `(f₀, err)`

of the estimated `f(x0)`

and an error estimate.

The return value of `f`

can be any type supporting `±`

and `norm`

operations (i.e. a normed vector space).
Similarly, `h`

and `x0`

can be in any normed vector space,
in which case `extrapolate`

performs Richardson extrapolation
of `f(x0+s*h)`

to `s=0⁺`

(i.e. it takes the limit as `x`

goes
to `x0`

along the `h`

direction).

On each step of Richardson extrapolation, it shrinks `x-x0`

by
a factor of `contract`

, stopping when the estimated error is
`< max(rtol*norm(f₀), atol)`

, when the estimated error
increases by more than `breaktol`

(e.g. due to numerical errors in the
computation of `f`

), when `f`

returns a non-finite value (`NaN`

or `Inf`

),
or when `f`

has been evaluated `maxeval`

times. Note that
if the function may converge to zero, you may want
specify a nonzero `atol`

(which cannot be set by default
because it depends on the scale/units of `f`

); alternatively,
in such cases `extrapolate`

will halt when it becomes
limited by the floating-point precision. (Passing `breaktol=Inf`

can be useful to force `extrapolate`

to continue shrinking `h`

even
if polynomial extrapolation is initially failing to converge,
possibly at the cost of extraneous function evaluations.)

If `x0 = ±∞`

(`±Inf`

), then `extrapolate`

computes the limit of
`f(x)`

as `x ⟶ ±∞`

using geometrically *increasing* values
of `h`

(by factors of `1/contract`

).

In general, the starting `h`

should be large enough that `f(x0+h)`

can be computed accurately and efficiently (e.g. without
severe cancellation errors), but small enough that `f`

does not
oscillate much between `x0`

and `x0+h`

. i.e. `h`

should be a typical
scale over which the function `f`

varies significantly.

Technically, Richardson extrapolation assumes that `f(x0+h)`

can
be expanded in a power series in `h^power`

, where the default
`power=1`

corresponds to an ordinary Taylor series (i.e. assuming
`f`

is analytic at `x0`

). If this is not true, you may obtain
slow convergence from `extrapolate`

, but you can pass a different
value of `power`

(e.g. `power=0.5`

) if your `f`

has some different
(Puiseux) power-series expansion. Conversely, if `f`

is
an *even* function around `x0`

, i.e. `f(x0+h) == f(x0-h)`

,
so that its Taylor series contains only *even* powers of `h`

,
you can accelerate convergence by passing `power=2`

.

`extrapolate(fh_itr; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf)`

Similar to `extrapolate(f, h)`

, performs Richardson extrapolation of a sequence of
values `f(h)`

to `h → 0`

, but takes an iterable collection `fh_itr`

of a
sequence of `(f(h), h)`

tuples (in order of decreasing `|h|`

).

There is no `contract`

keyword argument since the contraction factors are determined
by the sequence of `h`

values (which need not contract by the same amount). The
tolerances `atol`

and `rtol`

both default to `0`

so that by default it examines
*all* of the values in the `fh_itr`

collection. Otherwise, the keyword arguments
have the same meanings as in `extrapolate(f, h)`

.

`extrapolate!(fh::AbstractVector; power=1, atol=0, rtol=0, maxeval=typemax(Int), breaktol=Inf)`

Similar to `extrapolate(fh_itr)`

, performs Richardson extrapolation on an array `fh`

of `(f(h), h)`

tuples (in order of decreasing `|h|`

), but overwrites the array
`fh`

in-place with intermediate calculations.

(Thus, the array `fh`

must be a vector of `Tuple{T,H}`

values, where `H<:Number`

is
the type of `h`

and `T`

is the type of the extrapolated `f(0)`

**result**. This `T`

should be a floating-point type, i.e. `fh`

should contain `float(f(h))`

if the
function you are extrapolating is not already floating-point-valued.)

For example, let's extrapolate `sin(x)/x`

to `x=0`

(where the correct answer is `1`

) starting at `x=1`

, printing out the `x`

value at each step so that we can see what the algorithm is doing.

(Since `f`

is passed as the first argument to `extrapolate`

, we
can use Julia's `do`

syntax to conveniently define a multi-line
anonymous function to pass.)

```
extrapolate(1.0, rtol=1e-10) do x
@show x
sin(x)/x
end
```

giving the output:

```
x = 1.0
x = 0.125
x = 0.015625
x = 0.001953125
x = 0.000244140625
x = 3.0517578125e-5
(1.0000000000000002, 2.0838886172214188e-13)
```

That is, it evaluates our function `sin(x)/x`

for 6 different values of `x`

and returns `1.0000000000000002`

, which is accurate to machine precision (the error is `≈ 2.2e-16`

). The returned error estimate of `2e-13`

is conservative, which is typical for extrapolating well-behaved functions.

Since `sin(x)/x`

is an *even* (symmetric) function around `x=0`

,
its Taylor series contains only even powers of `x`

. We can
exploit this fact to *accelerate convergence for even functions* by
passing `power=2`

to `extrapolate`

:

```
extrapolate(1.0, rtol=1e-10, power=2) do x
@show x
sin(x)/x
end
```

gives

```
x = 1.0
x = 0.125
x = 0.015625
x = 0.001953125
x = 0.000244140625
(1.0, 0.0)
```

which converged to machine precision (in fact, the exact result) in only 5 function evaluations (1 fewer than above).

Using the `x0`

keyword argument, we can compute the limit of `f(x)`

as `x ⟶ x0`

. In fact, you can pass `x0 = Inf`

to compute a limit as
`x ⟶ ∞`

(which is accomplished internally by a change of variables `x = 1/u`

and performing Richardson extrapolation to `u=0`

). For example:

```
extrapolate(1.0, x0=Inf) do x
@show x
(x^2 + 3x - 2) / (x^2 + 5)
end
```

gives

```
x = 1.0
x = 8.0
x = 64.0
x = 512.0
x = 4096.0
x = 32768.0
x = 262144.0
(1.0000000000000002, 1.2938539128981574e-12)
```

which is the correct result (`1.0`

) to machine precision.

One nice application of infinite limits is extrapolating infinite series. If we start with an integer `x`

, then the default `contract=0.125`

will increase `x`

by a factor of `8.0`

on each iteration, so `x`

will always be an exact integer and we can use it as the number of terms in a series.

For example, suppose we are computing the infinite series `1/1² + 1/2² + 1/3² + ⋯`

. This is the famous Basel problem, and it converges to `π²/6 ≈ 1.644934066848226…`

. If we compute it by brute force, however, we need quite a few terms to get high accuracy:

```
julia> sum(n -> 1/n^2, 1:100) - π^2/6
-0.009950166663333482
julia> sum(n -> 1/n^2, 1:10^4) - π^2/6
-9.99950001654426e-5
julia> sum(n -> 1/n^2, 1:10^9) - π^2/6
-9.999985284281365e-10
```

Even with 10⁹ terms we get only about 9 digits. Instead, we can use `extrapolate`

(starting at 1 term):

```
julia> val, err = extrapolate(1, x0=Inf) do N
@show N
sum(n -> 1/n^2, 1:Int(N))
end
N = 1.0
N = 8.0
N = 64.0
N = 512.0
N = 4096.0
N = 32768.0
(1.6449340668482288, 4.384936858059518e-12)
julia> (val - π^2/6)/(π^2/6)
1.4848562646983628e-15
```

By `32768`

terms, the extrapolated value is accurate to about 15 digits.

A classic application of Richardson extrapolation is the accurate evaluation of derivatives via finite-difference approximations (although analytical derivatives, e.g. by automatic differentiation, are typically vastly more efficient when they are available). In this example, we use Richardson extrapolation on the forward-difference approximation `f'(x) ≈ (f(x+h)-f(x))/h`

, for which the error decreases as `O(h)`

but a naive application to a very small `h`

will yield a huge cancellation error from floating-point roundoff effects. We differentiate `f(x)=sin(x)`

at `x=1`

, for which the correct answer is `cos(1) ≈ 0.5403023058681397174009366...`

, starting with `h=0.1`

```
extrapolate(0.1, rtol=0) do h
@show h
(sin(1+h) - sin(1)) / h
end
```

Although we gave an `rtol`

of `0`

, the `extrapolate`

function will terminate after a finite number of steps when it detects that improvements are limited by floating-point error:

```
h = 0.1
h = 0.0125
h = 0.0015625
h = 0.0001953125
h = 2.44140625e-5
h = 3.0517578125e-6
(0.5403023058683176, 1.7075230118734908e-12)
```

The output `0.5403023058683176`

differs from `cos(1)`

by `≈ 1.779e-13`

, so in this case the returned error estimate is only a little conservative. Unlike the `sin(x)/x`

example, `extrapolate`

is not able
to attain machine precision (the floating-point cancellation error in this function is quite severe for small `h`

!), but it is able to get surprisingly close.

Another possibility for a finite-difference/Richardson combination was suggested by Ridders (1982), who computed both `f'(x)`

and `f''(x)`

(the first and second derivatives) simultaneously using a center-difference approximation, which requires two new `f(x)`

evaluations for each `h`

. In particular, the center-difference approximations are `f'(x) ≈ (f(x+h)-f(x-h))/2h`

and `f''(x) ≈ (f(x+h)-2f(x)+f(x-h))/h²`

, both of which have errors that go as `O(h²)`

. We can plug both of these functions *simultaneously* into `extrapolate`

(so that they share `f(x±h)`

evaluations) by using a vector-valued function returning `[f', f'']`

. Moreover, since these center-difference approximations are even functions of `h`

(identical for `±h`

), we can pass `power=2`

to `extrapolate`

in order to exploit the even-power Taylor expansion. Here is a function implementing both of these ideas:

```
# returns (f'(x), f''(x))
function riddersderiv2(f, x, h; atol=0, rtol=atol>0 ? 0 : sqrt(eps(typeof(float(real(x+h))))), contract=0.5)
f₀ = f(x)
val, err = extrapolate(h, atol=atol, rtol=rtol, contract=contract, power=2) do h
f₊, f₋ = f(x+h), f(x-h)
[(f₊-f₋)/2h, (f₊-2f₀+f₋)/h^2]
end
return val[1], val[2]
end
```

(This code could be made even more efficient by using StaticArrays.jl for the `[f', f'']`

vector.) The original paper by Ridders accomplishes something similar in `< 20`

lines of TI-59 calculator code, by the way; so much for high-level languages!

For example,

```
julia> riddersderiv2(1, 0.1, rtol=0) do x
@show x
sin(x)
end
x = 1
x = 1.1
x = 0.9
x = 1.05
x = 0.95
x = 1.025
x = 0.975
x = 1.0125
x = 0.9875
x = 1.00625
x = 0.99375
(0.5403023058681394, -0.841470984807975)
```

evaluates the first and second derivatives of `sin(x)`

at `x=1`

and obtains the correct answer `(cos(1), -sin(1))`

to about 15 and 13 decimal digits, respectively, using 11 function evaluations.

It is useful to consider a finite-difference approximation for the derivative of
the function `1/x`

at some `x ≠ 0`

: i.e. computing the limit of `f(h) = (1/(x+h) - 1/x) / h`

as `h`

goes to zero similar to above.

This function `f(h)`

has a pole at `h=-x`

, i.e. `f(-x)`

blows up. This means
that the Taylor series of `f(h)`

only converges for `h`

values small enough to avoid this pole, and in
fact for `|h| < |x|`

. Since Richardson extrapolation is essentially
approximating the Taylor series, this means that the extrapolation process doesn't
converge if the starting `h`

is too large, and `extrapolation`

will give up and
halt with the wrong answer.

This lack of convergence is easily observed: set `x=0.01`

(where the correct derivative of `1/x`

is `-10000`

) and consider what happens for a starting `h`

that is too large compared to `|x|`

:

```
julia> extrapolate(1.0) do h
@show h
x = 0.01
(1/(x+h) - 1/x) / h
end
h = 1.0
h = 0.125
h = 0.015625
(-832.4165749908325, 733.4066740007335)
```

Before reaching an `|h| < 0.01`

where the power series could begin to converge, `extrapolate`

gave up and returned a wrong answer (with a large error estimate to let you know that the result is garbage)! In contrast, if we start with a small enough `h`

then it converges just fine and returns the correct answer (`-10000`

) to nearly machine precision:

```
julia> extrapolate(0.01) do h
@show h
x = 0.01
(1/(x+h) - 1/x) / h
end
h = 0.01
h = 0.00125
h = 0.00015625
h = 1.953125e-5
h = 2.44140625e-6
h = 3.0517578125e-7
(-10000.000000000211, 4.066770998178981e-6)
```

Of course, if you know that your function blows up like this, it is easy to choose good initial `h`

, but how can we persuade `extrapolate`

to do a better job automatically?

The trick is to use the `breaktol`

keyword argument. `breaktol`

defaults to `2`

,
which means that `extrapolate`

gives up if the best error estimate increases by
more than a factor of 2 from one iteration to the next. Ordinarily, this
kind of breakdown in convergence arises because you hit the limits of floating-point precision, and halting extrapolation is the right thing to do. Here, however, it would converge if we just continued shrinking `h`

. So, we simply set `breaktol=Inf`

to force extrapolation to continue, which works even for a large initial `h=1000.0`

:

```
julia> extrapolate(1000.0, breaktol=Inf) do h
@show h
x = 0.01
(1/(x+h) - 1/x) / h
end
h = 1000.0
h = 125.0
h = 15.625
h = 1.953125
h = 0.244140625
h = 0.030517578125
h = 0.003814697265625
h = 0.000476837158203125
h = 5.9604644775390625e-5
h = 7.450580596923828e-6
h = 9.313225746154785e-7
h = 1.1641532182693481e-7
(-10000.000000029328, 5.8616933529265225e-8)
```

Not that it continues extrapolating until it reaches small `h`

values where the power-series converges, and in the end it again returns the correct answer to nearly machine precision (and would reach machine precision if we set a smaller `rtol`

). (The `extrapolate`

function *automatically* discards the initial points where the polynomial extrapolation fails.)