Numerical differentiation

Short description: Use of numerical analysis to estimate derivatives of functions

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. Unlike analytical differentiation, which provides exact expressions for derivatives, numerical differentiation relies on the function's values at a set of discrete points to estimate the derivative's value at those points or at intermediate points. This approach is particularly useful when dealing with data obtained from experiments, simulations, or situations where the function is defined only at specific intervals.

## Contents

* 1 Finite differences
* 2 Step size
* 2.1 Example
* 3 Three Point methods
* 3.1 Numerical Example
* 3.2 Code
* 4 Higher derivatives
* 4.1 Numerical Example
* 4.2 Arbitrary Derivatives
* 5 Complex-variable methods
* 6 Differential quadrature
* 7 See also
* 8 References
* 9 External links

## 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 (x, f(x)) and (x + h, f(x + h)).[1] Choosing a small number h, h represents a small change in x, and it can be either positive or negative. The slope of this line is f(x+h)−f(x)h. This expression is Newton's difference quotient (also known as a first-order divided difference).

To obtain an error estimate for this approximation, one can use Taylor expansion of f(x) about the base point x to give f(x+h)=f(x)+hf′(x)+h22f″(c) for some c between x and x+h. Rearranging gives f′(x)=f(x+h)−f(x)h⏟Slope of secant line−h2f″(c)⏟Error term. The slope of this secant line differs from the slope of the tangent line by an amount that is approximately proportional to h. As h approaches zero, the slope of the secant line approaches the slope of the tangent line and the error term vanishes. Therefore, the true derivative of f at x is the limit of the value of the difference quotient as the secant lines get closer and closer to being a tangent line: f′(x)=limh->0f(x+h)−f(x)h.

Since immediately substituting 0 for h results in 00 indeterminate form, calculating the derivative directly can be unintuitive.

Equivalently, the slope could be estimated by employing positions x − h and x.

Another two-point formula is to compute the slope of a nearby secant line through the points (x − h, f(x − h)) and (x + h, f(x + h)). The slope of this line is 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 h2. Hence for small values of h this is a more accurate approximation to the tangent line than the one-sided estimation. However, although the slope is being computed at x, the value of the function at x is not involved.

The estimation error is given by R=−f(3)(c)6h2, 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 h = 0.001.[2][3]

## Step size

See also: Adaptive step size

Example showing the difficulty of choosing h 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, h. To illustrate, consider the two-point approximation formula with error term:

f′(x)=f(x+h)−f(x−h)2h−f(3)(c)6h2

where c is some point between x−h and x+h. Let e(x) denote the roundoff error encountered when evaluating the function f(x) and f^(x) denote the computed value of f(x). Therefore, f(x)=f^(x)+e(x). The total error in the approximation is f′(x)−f^(x+h)−f^(x−h)2h=e(x+h)−e(x−h)2h⏟Roundoff error−f(3)(c)6h2⏟Truncation error. Assuming that the roundoff errors are bounded by some number ϵ>0 and the third derivative of f(x) is bounded by some number M>0, we get |f′(x)−f^(x+h)−f^(x−h)2h|≤ϵh+h26M.

To reduce the truncation error h26M, we must reduce h. But as h is reduced, the roundoff error ϵh increases. Due to the need to divide by the small value h, all the finite-difference formulae for numerical differentiation are similarly ill-conditioned.[4]

If instead the derivative is approximated using f′(x)=f(x+h)−f(x)h and we assume the relative approximation error due to rounding is bounded by machine epsilon and the values |f″(x)| and |f(x)| are bounded on [x0,x0+h] by M1and M2 respectively it can be shown that:[5]|f′(x0)−f^(x0+h)−f^(x0)h|≤h2M1+2ϵhM2By minimizing this upper bound an estimate for the optimal step size can be obtained. However, employing this method would require the knowledge of the bounds M1 and M2. Instead, if an approximate upper bound is used - one where M1and M2 are replaced by |f″(x0)| and |f(x0)| the error can be estimated as:e(h)=h2|f″(x0)|+2ϵh|f(x0)|By minimizing e(h) an estimate for the optimal step size h can be found:[5] h=2ε|f(x)f″(x)|(though not when f″(x)=0). If no information about the function f is known, sometimes the approximation |f(x)f″(x)|≈x is used giving a choice for h that is small without producing a large rounding error: εx, where the machine epsilon ε is typically of the order of 2.2×10−16 for double precision.[6] In practice, one cannot compute the value that minimizes the above error bounds and therefore cannot compute an approximation of the optimal step size without information about higher-order derivatives.[7]

### Example

To demonstrate this difficulty, consider approximating the derivative of the function f(x)=2x1+x at the point x0=9. In this case, we can calculate f′(x)=2+x(1+x)2 which gives f′(9)=2+9(1+9)2=2+3(1+3)2=516=0.3125. Using 64-bit floating point numbers, the following approximations are generated with the two-point approximation formula and increasingly smaller step sizes. The smallest absolute error is produced for a step size of 10−4, after which the absolute error steadily increases as the roundoff errors dominate calcuations.

Step Size (h) | Approximation | Absolute Error

|
|

10−1 | 0.31250397877426117 | 3.97877×10−6
10−2 | 0.31250003978589014 | 3.97858×10−8
10−3 | 0.3125000003976197 | 3.97619×10−10
10−4 | 0.31250000000593303 | 5.93303×10−12
10−5 | 0.3124999999215561 | 7.84439×10−11
10−6 | 0.31249999965510256 | 3.44897×10−10
10−7 | 0.31249999921101335 | 7.88986×10−10
10−8 | 0.31250002585636594 | 2.58563×10−8
10−9 | 0.31250024790097086 | 2.47900×10−7
10−10 | 0.31249669518729206 | 3.30481×10−6
10−11 | 0.31246116805050406 | 3.88319×10−5
10−12 | 0.312194714524594 | 3.05285×10−4
10−13 | 0.31530333899354446 | 2.80333×10−3
10−14 | 3.552713678800501 | 4.27713×10−2

For computer calculations the problems are exacerbated because, although x necessarily holds a representable floating-point number in some precision (32 or 64-bit, etc.), x \+ h almost certainly will not be exactly representable in that precision. This means that x \+ h will be changed (by rounding or truncation) to a nearby machine-representable number, with the consequence that (x + h) − x will not equal h; the two function evaluations will not be exactly h 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 h = 0.1 will not be a round number in binary; it is 0.000110011001100...2 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 dx and h are the same. With C and similar languages, a directive that xph is a volatile variable will prevent this.

## Three Point methods

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

f(x+h)=f(x)+hf′(x)+h22!f″(x)+h33!f‴(x)+...

|
|

(1)

Replacing h by 2h gives

f(x+2h)=f(x)+2hf′(x)+4h22!f″(x)+8h33!f‴(x)+...

|
|

(2)

Multiplying identity (1) by 4 gives

4f(x+h)=4f(x)+4hf′(x)+4h22!f″(x)+4h33!f‴(x)+...

|
|

(1')

Subtracting identity (1') from (2) eliminates the h2 term:

f(x+2h)−4f(x+h)=−3f(x)−2hf′(x)+4h33!f‴(x)+...

which can be written as

f(x+2h)−4f(x+h)=−3f(x)−2hf′(x)+O(h3).

Rearranging terms gives f′(x)=−3f(x)+4f(x+h)−f(x+2h)2h+O(h2),

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

f′(x)=f(x+h)−f(x−h)2h+O(h2) which is called the three-point central difference formula, and f′(x)=f(x−2h)−4f(x−h)+3f(x)2h+O(h2) which is called the three-point backward difference formula.

By a similar approach, the five point midpoint approximation formula can be derived as:[8] f′(x)=−f(x+2h)+8f(x+h)−8f(x−h)+f(x−2h)12h+O(h4).

### Numerical Example

Consider approximating the derivative of f(x)=xsin⁡x at the point x0=π4. Since f′(x)=sin⁡x+xcos⁡x, the exact value is f′(π4)=sin⁡π4+π4cos⁡π4=12+π42≈1.2624671484563432.

Formula  | Step Size (h)  | Approximation  | Absolute Error

|
|
|

Three-point forward difference formula | 10−1 | 1.2719084899816118 | 9.441×10−3
10−2 | 1.2625569346253918 | 8.978×10−5
10−3 | 1.2624680412510747 | 8.927×10−7
10−4 | 1.2624671573796542 | 8.923×10−9
Three-point backward difference formula | 10−1 | 1.2719084899816118 | 8.307×10−3
10−2 | 1.2625569346253918 | 8.864×10−5
10−3 | 1.2624680412510747 | 8.916×10−7
10−4 | 1.2624671573796542 | 8.922×10−9
Three-point central difference formula | 10−1 | 1.2719084899816118 | 4.457×10−3
10−2 | 1.2625569346253918 | 4.461×10−5
10−3 | 1.2624680412510747 | 4.461×10−7
10−4 | 1.2624671573796542 | 4.461×10−9

### Code

The following is an example of a Python implementation for finding derivatives numerically for f(x)=2x1+x using the various three-point difference formulas at x0=4. The function `func` has derivative `func_prime`.

Example implementation in Python

import math

def func(x):
return (2*x) / (1 + math.sqrt(x))

def func_prime(x):
return (2 + math.sqrt(x)) / ((1 + math.sqrt(x))**2)

def three_point_forward(value, h):
return ((-3/2) * func(value) + 2*func(value + h) - (1/2)*func(value + 2*h)) / h

def three_point_central(value, h):
return ((-1/2)*func(value - h) + (1/2)*func(value + h)) / h

def three_point_backward(value, h):
return ((1/2)*func(value - 2*h) - 2*func(value - h) + (3/2)*func(value)) / h

value = 4
actual = func_prime(value)

print("Actual value " + str(actual))
print("============================================")

for step_size in [0.1, 0.01, 0.001, 0.0001]:

print("Step size " + str(step_size))

forward = three_point_forward(value, step_size)
backward = three_point_backward(value, step_size)
central = three_point_central(value, step_size)

print("Forward  {:>12}, Error = {:>12}".format(str(forward), str(abs(forward - actual))))
print("Backward {:>12}, Error = {:>12}".format(str(forward), str(abs(backward - actual))))
print("Central  {:>12}, Error = {:>12}".format(str(forward), str(abs(central - actual))))
print("============================================")

Output

Actual value 0.4444444444444444
============================================
Step size 0.1
Forward  0.4443963018050967, Error = 4.814263934771468e-05
Backward 0.4443963018050967, Error = 2.5082646145202503e-05
Central  0.4443963018050967, Error = 5.231976394060034e-05
============================================
Step size 0.01
Forward  0.4444439449793336, Error = 4.994651108258807e-07
Backward 0.4444439449793336, Error = 2.507721614808389e-07
Central  0.4444439449793336, Error = 5.036366184096863e-07
============================================
Step size 0.001
Forward  0.4444444394311464, Error = 5.013297998957e-09
Backward 0.4444444394311464, Error = 2.507574814458735e-09
Central  0.4444444394311464, Error = 5.017960935660426e-09
============================================
Step size 0.0001
Forward  0.4444444443896245, Error = 5.4819926376126205e-11
Backward 0.4444444443896245, Error = 2.5116131396885066e-11
Central  0.4444444443896245, Error = 5.037903427762558e-11
============================================

## Higher derivatives

Using Taylor Series, it is possible to derive formulas to approximate the second (and higher order) derivatives of a general function. For a function f(x) and some number h>0, expanding it about x+h and x−h gives

f(x+h)=f(x)+f′(x)h+12f″(x)h2+16f‴(x)h3+124f(4)(ξ1)h4 and f(x−h)=f(x)−f′(x)h+12f″(x)h2−16f‴(x)h3+124f(4)(ξ2)h4 where x−h<ξ2<x<ξ1<x+h. Adding these two equations gives

f(x+h)+f(x−h)=2f(x)+f″(x)h2+124[f(4)(ξ1)+f(4)(ξ2)]h4.

If f(4) is continuous on [x−h,x+h], then 12[f(4)(ξ1)+f(4)(ξ2)] is between f(4)(ξ1) and f(4)(ξ2). The Intermediate Value Theorem guarantees a number, say ξ, between ξ1 and ξ2. We can therefore write

f″(x)=1h2[f(x−h)−2f(x)+f(x+h)]−h212f(4)(ξ).

where x−h<ξ<x+h.

### Numerical Example

Consider approximating the second derivative of the function f(x)=exsin⁡x at the point x0=π4.

Since f″(x)=2excos⁡x, the exact value is f″(π4)=2eπ4≈3.1017663.

Step Size (h)  | Approximation  | Absolute Error

|
|

10−1 | 3.096593338003672 | 5.17305×10−3
10−2 | 3.1017146973844056 | 5.16964×10−5
10−3 | 3.1017658768117684 | 5.17024×10−7
10−4 | 3.1017663770782633 | 1.67577×10−8

### Arbitrary Derivatives

Using Newton's difference quotient, f′(x)=limh->0f(x+h)−f(x)h the following can be shown[9] (for n > 0): f(n)(x)=limh->01hn∑k=0n(−1)k+n(nk)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,[10] the first derivative can be calculated by the complex-step derivative formula:[11][12][13] f′(x)=ℑ(f(x+ih))h+O(h2),i2:=−1.

The recommended step size to obtain accurate derivatives for a range of conditions is h=10−200.[14] This formula can be obtained by Taylor series expansion: f(x+ih)=f(x)+ihf′(x)−12!h2f″(x)−i3!h3f(3)(x)+⋯.

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.[15][16][17] f(n)(x)≈𝒞n2−1(n)(f(x+i(1)h+⋯+i(n)h))hn where the i(k) denote the multicomplex imaginary units; i(1)≡i. The 𝒞k(n) operator extracts the kth component of a multicomplex number of level n, e.g., 𝒞0(n) extracts the real component and 𝒞n2−1(n) extracts the last, “most imaginary” component. The method can be applied to mixed derivatives, e.g. for a second-order derivative ∂2f(x,y)∂x∂y≈𝒞3(2)(f(x+i(1)h,y+i(2)h))h2

A C++ implementation of multicomplex arithmetics is available.[18]

In general, derivatives of any order can be calculated using Cauchy's integral formula:[19] f(n)(a)=n!2πi∮γf(z)(z−a)n+1dz, where the integration is done numerically.

Using complex variables for numerical differentiation was started by Lyness and Moler in 1967.[20] 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.[21] An algorithm that can be used without requiring knowledge about the method or the character of the function was developed by Fornberg.[4]

## Differential quadrature

Differential quadrature is the approximation of derivatives by using weighted sums of function values.[22][23] 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.[24]

## See also

* Automatic differentiation - Numerical calculations carrying along derivatives
* Five-point stencil - A point and its four nearest neighbors
* List of numerical-analysis software
* Numerical integration - Methods of calculating definite integrals
* Numerical methods for ordinary differential equations - Methods used to find numerical solutions of ordinary differential equations
* Savitzky–Golay filter - Algorithm to smooth data points

## References

1. ↑ Richard L. Burden, J. Douglas Faires (2000), Numerical Analysis, (7th Ed), Brooks/Cole. ISBN 0-534-38216-9.
2. ↑ Katherine Klippert Merseth (2003). Windows on Teaching Math: Cases of Middle and Secondary Classrooms. Teachers College Press. p. 34. ISBN 978-0-8077-4279-2. https://archive.org/details/windowsonteachin00mers.
3. ↑ Tamara Lefcourt Ruby; James Sellers; Lisa Korf; Jeremy Van Horn; Mike Munn (2014). Kaplan AP Calculus AB & BC 2015. Kaplan Publishing. p. 299. ISBN 978-1-61865-686-5.
4. ↑ 4.0 4.1 Numerical Differentiation of Analytic Functions, B Fornberg – ACM Transactions on Mathematical Software (TOMS), 1981.
5. ↑ 5.0 5.1 Mørken, Knut (September 2017). Numerical Algorithms and Digital Representation lecture notes. University of Oslo. pp. 269-290. https://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h21/kompendiet/matinf1100.pdf.
6. ↑ Following Numerical Recipes in C, Chapter 5.7.
7. ↑ Richard L. Burden, J. Douglas Faires (2000), Numerical Analysis, (7th Ed), Brooks/Cole. ISBN 0-534-38216-9.
8. ↑ Abramowitz & Stegun, Table 25.2.
9. ↑ Shilov, George. Elementary Real and Complex Analysis.
10. ↑ Using Complex Variables to Estimate Derivatives of Real Functions, W. Squire, G. Trapp – SIAM REVIEW, 1998.
11. ↑ Martins, J. R. R. A.; Sturdza, P.; Alonso, J. J. (2003). "The Complex-Step Derivative Approximation". ACM Transactions on Mathematical Software 29 (3): 245–262. doi:10.1145/838250.838251.
12. ↑ Differentiation With(out) a Difference by Nicholas Higham
13. ↑ article from MathWorks blog, posted by Cleve Moler
14. ↑ Martins, Joaquim R. R. A.; Ning, Andrew (2021-10-01) (in en). Engineering Design Optimization. Cambridge University Press. ISBN 978-1108833417. http://flowlab.groups.et.byu.net/mdobook.pdf.
15. ↑ "Archived copy". http://russell.ae.utexas.edu/FinalPublications/ConferencePapers/2010Feb_SanDiego_AAS-10-218_mulicomplex.pdf.
16. ↑ Lantoine, G.; Russell, R. P.; Dargent, Th. (2012). "Using multicomplex variables for automatic computation of high-order derivatives". ACM Trans. Math. Softw. 38 (3): 1–21. doi:10.1145/2168773.2168774.
17. ↑ Verheyleweghen, A. (2014). "Computation of higher-order derivatives using the multi-complex step method". http://folk.ntnu.no/preisig/HAP_Specials/AdvancedSimulation_files/2014/AdvSim-2014__Verheule_Adrian_Complex_differenetiation.pdf.
18. ↑ Bell, I. H. (2019). "mcx (multicomplex algebra library)". https://github.com/ianhbell/mcx.
19. ↑ Ablowitz, M. J., Fokas, A. S.,(2003). Complex variables: introduction and applications. Cambridge University Press. Check theorem 2.6.2
20. ↑ Lyness, J. N.; Moler, C. B. (1967). "Numerical differentiation of analytic functions". SIAM J. Numer. Anal. 4 (2): 202–210. doi:10.1137/0704019. Bibcode: 1967SJNA....4..202L.
21. ↑ Abate, J; Dubner, H (March 1968). "A New Method for Generating Power Series Expansions of Functions". SIAM J. Numer. Anal. 5 (1): 102–112. doi:10.1137/0705008. Bibcode: 1968SJNA....5..102A.
22. ↑ Differential Quadrature and Its Application in Engineering: Engineering Applications, Chang Shu, Springer, 2000, ISBN 978-1-85233-209-9.
23. ↑ Advanced Differential Quadrature Methods, Yingyan Zhang, CRC Press, 2009, ISBN 978-1-4200-8248-7.
24. ↑ Ahnert, Karsten; Abel, Markus (2007). "Numerical differentiation of experimental data: local versus global methods". Computer Physics Communications 177 (10): 764–774. doi:10.1016/j.cpc.2007.03.009. ISSN 0010-4655. Bibcode: 2007CoPhC.177..764A.

## External links

| Wikibooks has a book on the topic of: Numerical Methods

|

* 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.

0.00

(0 votes)

Original source: https://en.wikipedia.org/wiki/Numerical differentiation. Read more

|

Retrieved from "https://handwiki.org/wiki/index.php?title=Numerical_differentiation&oldid=4226689"