Numerical differentiation

Finite difference estimation of derivative

In numerical analysis, numerical differentiation algorithms estimate the derivative of a mathematical function or subroutine using values of the function and perhaps other knowledge about the function.

Finite differences

The simplest method is to use finite difference approximations.

A simple two-point estimation is to compute the slope of a nearby secant line through the points Template:Math and Template:Math. Choosing a small number Template:Mvar, Template:Mvar represents a small change in Template:Mvar, and it can be either positive or negative. The slope of this line is $$\frac{f(x+h)-f(x)}{h}.$$ This expression is Newton's difference quotient (also known as a first-order divided difference).

The slope of this secant line differs from the slope of the tangent line by an amount that is approximately proportional to Template:Mvar. As Template:Mvar approaches zero, the slope of the secant line approaches the slope of the tangent line. Therefore, the true derivative of Template:Math at Template:Mvar is the limit of the value of the difference quotient as the secant lines get closer and closer to being a tangent line: $$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h}.$$

Since immediately substituting 0 for Template:Mvar results in $\frac{0}{0}$ indeterminate form, calculating the derivative directly can be unintuitive.

Equivalently, the slope could be estimated by employing positions Template:Math and Template:Mvar.

Another two-point formula is to compute the slope of a nearby secant line through the points Template:Math and Template:Math. The slope of this line is $$\frac{f(x+h)-f(x-h)}{2h}.$$

This formula is known as the symmetric difference quotient. In this case the first-order errors cancel, so the slope of these secant lines differ from the slope of the tangent line by an amount that is approximately proportional to $h^2$. Hence for small values of Template:Mvar this is a more accurate approximation to the tangent line than the one-sided estimation. However, although the slope is being computed at Template:Mvar, the value of the function at Template:Mvar is not involved.

The estimation error is given by $$R=\frac{-f^{(3)}(c)}{6}{h}^2,$$ where $c$ is some point between $x-h$ and $x+h$. This error does not include the rounding error due to numbers being represented and calculations being performed in limited precision.

The symmetric difference quotient is employed as the method of approximating the derivative in a number of calculators, including TI-82, TI-83, TI-84, TI-85, all of which use this method with Template:Math.

Step size

Example showing the difficulty of choosing Template:Mvar due to both rounding error and formula error

An important consideration in practice when the function is calculated using floating-point arithmetic of finite precision is the choice of step size, Template:Mvar. If chosen too small, the subtraction will yield a large rounding error. In fact, all the finite-difference formulae are ill-conditioned and due to cancellation will produce a value of zero if Template:Mvar is small enough. If too large, the calculation of the slope of the secant line will be more accurately calculated, but the estimate of the slope of the tangent by using the secant could be worse.

For basic central differences, the optimal step is the cube-root of machine epsilon. For the numerical derivative formula evaluated at Template:Mvar and Template:Math, a choice for Template:Mvar that is small without producing a large rounding error is $\sqrt{\epsilon }x$ (though not when x = 0), where the machine epsilon Template:Mvar is typically of the order of Template:Val for double precision. A formula for Template:Mvar that balances the rounding error against the secant error for optimum accuracy is $$h=2\sqrt{\epsilon \left|\frac{f(x)}{f^{″}(x)}\right|}$$ (though not when $f^{″}(x)=0$), and to employ it will require knowledge of the function.

For computer calculations the problems are exacerbated because, although Template:Mvar necessarily holds a representable floating-point number in some precision (32 or 64-bit, etc.), Template:Math almost certainly will not be exactly representable in that precision. This means that Template:Math will be changed (by rounding or truncation) to a nearby machine-representable number, with the consequence that Template:Math will not equal Template:Mvar; the two function evaluations will not be exactly Template:Mvar apart. In this regard, since most decimal fractions are recurring sequences in binary (just as 1/3 is in decimal) a seemingly round step such as Template:Math will not be a round number in binary; it is 0.000110011001100...₂ A possible approach is as follows:

h     := sqrt(eps) * x;
xph   := x + h;
dx    := xph - x;
slope := (F(xph) - F(x)) / dx;

However, with computers, compiler optimization facilities may fail to attend to the details of actual computer arithmetic and instead apply the axioms of mathematics to deduce that Template:Math and Template:Mvar are the same. With C and similar languages, a directive that Template:Math is a volatile variable will prevent this.

Other methods

Higher-order methods

Template:Further To obtain more general derivative approximation formulas for some function $f(x)$, let $h>0$ be a positive number close to zero. The Taylor expansion of $f(x)$ about the base point $x$ is Template:NumBlk

Replacing $h$ by $2h$ gives Template:NumBlk

Multiplying identity (Template:EquationNote) by 4 gives Template:NumBlk

Subtracting identity (Template:EquationNote) from (Template:EquationNote) eliminates the $h^2$ term:

$$f(x+2h)-4f(x+h)=-3f(x)-2h{f}^{\prime }(x)+\frac{4h^3}{3!}{f}^{‴}(x)+...$$

which can be written as

$$f(x+2h)-4f(x+h)=-3f(x)-2h{f}^{\prime }(x)+O(h^3).$$

Rearranging terms gives $$f^{\prime }(x)=\frac{-3f(x)+4f(x+h)-f(x+2h)}{2h}+O(h^2),$$

which is called the three-point forward difference formula for the derivative. Using a similar approach, one can show

$$f^{\prime }(x)=\frac{f(x+h)-f(x-h)}{2h}+O(h^2)$$ which is called the three-point central difference formula, and $$f^{\prime }(x)=\frac{f(x-2h)-4f(x-h)+3f(x)}{2h}+O(h^2)$$ which is called the three-point backward difference formula.

By a similar approach, the five point midpoint approximation formula can be derived as: $$f^{\prime }(x)=\frac{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}{12h}+O(h^4).$$

Higher derivatives

Using Newton's difference quotient, $$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h}$$ the following can be shown (for Template:Math): $$f^{(n)}(x)=\underset{h\to 0}{\lim }\frac{1}{h^n}\sum _{k=0}^n(-1{)}^{k+n}(\genfrac{}{}{0ex}{}{n}{k})f(x+kh)$$

Complex-variable methods

The classical finite-difference approximations for numerical differentiation are ill-conditioned. However, if $f$ is a holomorphic function, real-valued on the real line, which can be evaluated at points in the complex plane near $x$, then there are stable methods. For example, the first derivative can be calculated by the complex-step derivative formula: $$f^{\prime }(x)=\frac{\mathrm{\Im }(f(x+\mathrm{i}h))}{h}+O(h^2),{\mathrm{i}}^{\mathrm{2}}:=-1.$$

The recommended step size to obtain accurate derivatives for a range of conditions is $h=10^{-200}$. This formula can be obtained by Taylor series expansion: $$f(x+\mathrm{i}h)=f(x)+\mathrm{i}h{f}^{\prime }(x)-\frac{1}{2!}{h}^2{f}^{″}(x)-\frac{\mathrm{i}}{3!}{h}^3{f}^{(3)}(x)+\cdots .$$

The complex-step derivative formula is only valid for calculating first-order derivatives. A generalization of the above for calculating derivatives of any order employs multicomplex numbers, resulting in multicomplex derivatives. $$f^{(n)}(x)\approx \frac{{𝒞}_{n^2-1}^{(n)}(f(x+{\mathrm{i}}^{(1)}h+\cdots +{\mathrm{i}}^{(n)}h))}{h^n}$$ where the ${\mathrm{i}}^{(k)}$ denote the multicomplex imaginary units; ${\mathrm{i}}^{(1)}\equiv \mathrm{i}$. The ${𝒞}_k^{(n)}$ operator extracts the $k$th component of a multicomplex number of level $n$, e.g., ${𝒞}_0^{(n)}$ extracts the real component and ${𝒞}_{n^2-1}^{(n)}$ extracts the last, “most imaginary” component. The method can be applied to mixed derivatives, e.g. for a second-order derivative $$\frac{{\partial }^{2}f(x,y)}{\partial x\partial y}\approx \frac{{𝒞}_3^{(2)}(f(x+{\mathrm{i}}^{(1)}h,y+{\mathrm{i}}^{(2)}h))}{h^2}$$

A C++ implementation of multicomplex arithmetics is available.

In general, derivatives of any order can be calculated using Cauchy's integral formula: $$f^{(n)}(a)=\frac{n!}{2\pi i}{\oint }_{\gamma }\frac{f(z)}{(z-a{)}^{n+1}}\mathrm{d}z,$$ where the integration is done numerically.

Using complex variables for numerical differentiation was started by Lyness and Moler in 1967. Their algorithm is applicable to higher-order derivatives.

A method based on numerical inversion of a complex Laplace transform was developed by Abate and Dubner. An algorithm that can be used without requiring knowledge about the method or the character of the function was developed by Fornberg.

Differential quadrature

Differential quadrature is the approximation of derivatives by using weighted sums of function values. Differential quadrature is of practical interest because its allows one to compute derivatives from noisy data. The name is in analogy with quadrature, meaning numerical integration, where weighted sums are used in methods such as Simpson's rule or the trapezoidal rule. There are various methods for determining the weight coefficients, for example, the Savitzky–Golay filter. Differential quadrature is used to solve partial differential equations. There are further methods for computing derivatives from noisy data.

See also

• Template:Annotated link
• Template:Annotated link
• List of numerical-analysis software
• Template:Annotated link
• Template:Annotated link
• Template:Annotated link

External links

• Numerical Differentiation from wolfram.com
• NAG Library numerical differentiation routines
• Boost. Math numerical differentiation, including finite differencing and the complex step derivative
• Differentiation With(out) a Difference by Nicholas Higham, SIAM News.