View source for Numerical differentiation

{{Short description|Use of numerical analysis to estimate derivatives of functions}}
[[File:Derivative.svg|thumb|[[Finite difference]] estimation of derivative|230px|right]]
In [[numerical analysis]], '''numerical differentiation''' [[algorithm]]s estimate the [[derivative]] of a [[Function (mathematics)|mathematical function]] or [[Function (computer programming)|subroutine]] using values of the function and perhaps other knowledge about the function.

==Finite differences==
{{further|Finite difference}}

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 {{math|(''x'', ''f''(''x''))}} and {{math|(''x'' + ''h'', ''f''(''x'' + ''h''))}}.<ref>Richard L. Burden, J. Douglas Faires (2000), ''Numerical Analysis'', (7th Ed),  Brooks/Cole. {{isbn|0-534-38216-9}}.</ref> Choosing a small number {{mvar|h}}, {{mvar|h}} represents a small change in {{mvar|x}}, and it can be either positive or negative.  The slope of this line is
<math display="block">\frac{f(x + h) - f(x)}{h}.</math>
This expression is [[Isaac Newton|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 {{mvar|h}}. As {{mvar|h}} approaches zero, the slope of the secant line approaches the slope of the tangent line. Therefore, the true '''derivative of {{math|''f''}} at {{mvar|x}}''' is the limit of the value of the difference quotient as the secant lines get closer and closer to being a tangent line:
<math display="block">f'(x) = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h}.</math>

Since immediately [[substitution (logic)|substituting]] 0 for {{mvar|h}} results in <math>\frac{0}{0}</math> [[indeterminate form]], calculating the derivative directly can be unintuitive.

Equivalently, the slope could be estimated by employing positions {{math|''x'' − ''h''}} and {{mvar|x}}.

Another two-point formula is to compute the slope of a nearby secant line through the points {{math|(''x'' − ''h'', ''f''(''x'' − ''h''))}} and {{math|(''x'' + ''h'', ''f''(''x'' + ''h''))}}. The slope of this line is
<math display="block">\frac{f(x + h) - f(x - h)}{2h}.</math>

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 <math>h^2</math>. Hence for small values of {{mvar|h}} this is a more accurate approximation to the tangent line than the one-sided estimation. However, although the slope is being computed at {{mvar|x}}, the value of the function at {{mvar|x}} is not involved.

The estimation error is given by
<math display="block">R = \frac{-f^{(3)}(c)}{6} h^2,</math>
where <math>c</math> is some point between <math>x - h</math> and <math>x + h</math>.
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 {{math|1=''h'' = 0.001}}.<ref name="Merseth2003">{{cite book |author=Katherine Klippert Merseth |title=Windows on Teaching Math: Cases of Middle and Secondary Classrooms |url=https://archive.org/details/windowsonteachin00mers |url-access=limited |year=2003 |publisher=Teachers College Press |isbn=978-0-8077-4279-2 |page=[https://archive.org/details/windowsonteachin00mers/page/n60 34]}}</ref><ref name="RubySellers2014">{{cite book |author1=Tamara Lefcourt Ruby |author2=James Sellers |author3=Lisa Korf |author4=Jeremy Van Horn |author5=Mike Munn |title=Kaplan AP Calculus AB & BC 2015 |year=2014 |publisher=Kaplan Publishing |isbn=978-1-61865-686-5 |page=299}}</ref>

==Step size==
{{see also|Adaptive step size}}
[[Image:AbsoluteErrorNumericalDifferentiationExample.png|thumb|300px|Example showing the difficulty of choosing {{mvar|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, {{mvar|h}}. If chosen too small, the subtraction will yield a large [[rounding error]]. In fact, all the finite-difference formulae are [[ill-conditioned]]<ref name=Fornberg1>Numerical Differentiation of Analytic Functions, B Fornberg – ACM Transactions on Mathematical Software (TOMS), 1981.</ref> and due to cancellation will produce a value of zero if {{mvar|h}} is small enough.<ref name=SquireTrapp1>Using Complex Variables to Estimate Derivatives of Real Functions, W. Squire, G. Trapp – SIAM REVIEW, 1998.</ref> 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.<ref name="edo2021">{{Cite book|url= http://flowlab.groups.et.byu.net/mdobook.pdf | title=Engineering Design Optimization | last1=Martins|first1=Joaquim R. R. A. | last2=Ning|first2=Andrew | date=2021-10-01 | publisher=Cambridge University Press | isbn=978-1108833417|language=en}}</ref>

For basic central differences, the optimal step is the cube-root of [[machine epsilon]].<ref>Sauer, Timothy (2012). ''Numerical Analysis''. Pearson. p.248.</ref>
For the numerical derivative formula evaluated at {{mvar|x}} and {{math|''x'' + ''h''}}, a choice for {{mvar|h}} that is small without producing a large rounding error is <math>\sqrt{\varepsilon} x</math> (though not when ''x'' = 0), where the [[machine epsilon]] {{mvar|ε}} is typically of the order of {{val|2.2e−16}} for [[double precision]].<ref>Following ''[[Numerical Recipes]] in C'', [http://www.nrbook.com/a/bookcpdf/c5-7.pdf Chapter 5.7].</ref> A formula for {{mvar|h}} that balances the rounding error against the secant error for optimum accuracy is<ref>[http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h10/kompendiet/kap11.pdf p. 263].</ref>
<math display="block">h = 2\sqrt{\varepsilon\left|\frac{f(x)}{f''(x)}\right|}</math>
(though not when <math>f''(x) = 0</math>), and to employ it will require knowledge of the function.

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

==Other methods==
===Higher-order methods===
{{further|Finite difference coefficient}}
To obtain more general derivative approximation formulas for some function <math>f(x)</math>, let <math>h>0</math> be a positive number close to zero. The Taylor expansion of <math>f(x)</math> about the base point <math>x</math> is
{{NumBlk|::|<math>f(x+h) = f(x) + hf'(x) + \frac{h^{2}}{2!}f''(x) + \frac{h^{3}}{3!}f'''(x) + ...</math>|{{EquationRef|1}}}}

Replacing <math>h</math> by <math>2h</math> gives
{{NumBlk|::|<math>f(x+2h) = f(x) + 2hf'(x) + \frac{4h^{2}}{2!}f''(x) + \frac{8h^{3}}{3!}f'''(x) + ...</math>|{{EquationRef|2}}}}

Multiplying identity ({{EquationNote|1|}}) by 4 gives
{{NumBlk|::|<math>4f(x+h) = 4f(x) + 4hf'(x) + \frac{4h^{2}}{2!}f''(x) + \frac{4h^{3}}{3!}f'''(x) + ...</math>|{{EquationRef|1'}}}}

Subtracting identity ({{EquationNote|1'|}}) from ({{EquationNote|2|}}) eliminates the <math>h^{2}</math> term:

<math display="block">
f(x+2h) - 4f(x+h) = -3f(x) -2hf'(x) + \frac{4h^{3}}{3!}f'''(x) + ...</math>

which can be written as

<math display="block">
f(x+2h) - 4f(x+h) = -3f(x) -2hf'(x) + O(h^{3}).</math>

Rearranging terms gives
<math display="block">
f'(x) = \frac{-3f(x) + 4f(x+h) -f(x+2h)}{2h} + O(h^{2}),</math>

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

<math display="block">
f'(x) = \frac{f(x+h) - f(x-h)}{2h} + O(h^{2})</math>
which is called the '''three-point central difference formula''', and
<math display="block">
f'(x) = \frac{f(x-2h) -4f(x-h) + 3f(x)}{2h} + O(h^{2})</math>
which is called the '''three-point backward difference formula'''.

By a similar approach, the five point midpoint approximation formula can be derived as:<ref>Abramowitz & Stegun, Table 25.2.</ref>
<math display="block">f'(x) = \frac{-f(x + 2h) + 8 f(x + h) - 8 f(x - h) + f(x - 2h)}{12h} + O(h^{4}).</math>

=== Higher derivatives ===
Using Newton's difference quotient,
<math display="block">f'(x) = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h}</math>
the following can be shown<ref>{{cite book |last1=Shilov |first1=George |title=Elementary Real and Complex Analysis}}</ref> (for {{math|''n'' > 0}}):
<math display="block">f^{(n)}(x) = \lim_{h\to 0} \frac{1}{h^n} \sum_{k=0}^n (-1)^{k+n} \binom{n}{k} f(x + kh)</math>
<!-- [the following clearly has a wrong sign for n = 1]
However, a factor seems to be missing:
<math display="block">f^{(n)}(x)= (-1)^n \lim_{h\to 0} \frac{1}{h^n} \sum_{k=0}^n (-1)^{k+1} \binom{n}{k} f(x + kh)</math>
-->

==Complex-variable methods==

The classical finite-difference approximations for numerical differentiation are ill-conditioned. However, if <math>f</math> is a [[holomorphic function]], real-valued on the real line, which can be evaluated at points in the complex plane near <math>x</math>, then there are [[Numerical stability|stable]] methods. For example,<ref name=SquireTrapp1/> the first derivative can be calculated by the complex-step derivative formula:<ref>{{cite journal | last1 = Martins | first1 = J. R. R. A. | first2 = P. | last2 = Sturdza | first3 = J. J. | last3 = Alonso | year = 2003 | citeseerx=10.1.1.141.8002 | title = The Complex-Step Derivative Approximation | journal = ACM Transactions on Mathematical Software | volume = 29 | issue = 3 | pages = 245–262 | doi=10.1145/838250.838251| s2cid = 7022422 }}</ref><ref>[https://sinews.siam.org/Details-Page/differentiation-without-a-difference Differentiation With(out) a Difference] by [[Nicholas Higham]] </ref><ref>[https://blogs.mathworks.com/cleve/2013/10/14/complex-step-differentiation/ article] from [[MathWorks]] blog, posted by [[Cleve Moler]]</ref>
<math display="block">f'(x) = \frac{\Im(f(x + \mathrm{i}h))}{h} + O(h^2), \quad \mathrm{i^2}:=-1.</math>

The recommended step size to obtain accurate derivatives for a range of conditions is <math>h = 10^{-200}</math>.<ref name="edo2021"/>
This formula can be obtained by [[Taylor series]] expansion:
<math display="block">f(x+\mathrm{i}h) = f(x) + \mathrm{i}h f'(x) - \tfrac{1}{2!} h^2 f''(x) - \tfrac{\mathrm{i}}{3!} h^3 f^{(3)}(x) + \cdots.</math>

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.<ref>{{Cite web |url=http://russell.ae.utexas.edu/FinalPublications/ConferencePapers/2010Feb_SanDiego_AAS-10-218_mulicomplex.pdf |title=Archived copy |access-date=2012-11-24 |archive-url=https://web.archive.org/web/20140109145840/http://russell.ae.utexas.edu/FinalPublications/ConferencePapers/2010Feb_SanDiego_AAS-10-218_mulicomplex.pdf |archive-date=2014-01-09 |url-status=dead }}</ref><ref>{{cite journal | last1 = Lantoine | first1 = G. | last2 = Russell | first2 = R. P. | last3 = Dargent | first3 = Th. | title = Using multicomplex variables for automatic computation of high-order derivatives | journal = ACM Trans. Math. Softw. | volume = 38 | year = 2012 | issue = 3 | pages = 1–21 | doi = 10.1145/2168773.2168774 | s2cid = 16253562 }}</ref><ref>{{cite web | last1 = Verheyleweghen | first1 = A. | title = Computation of higher-order derivatives using the multi-complex step method | year = 2014 | url = http://folk.ntnu.no/preisig/HAP_Specials/AdvancedSimulation_files/2014/AdvSim-2014__Verheule_Adrian_Complex_differenetiation.pdf}}</ref>
<math display="block">f^{(n)}(x) \approx \frac{\mathcal{C}^{(n)}_{n^2-1}(f(x + \mathrm{i}^{(1)} h + \cdots + \mathrm{i}^{(n)} h))}{h^n}</math>
where the <math>\mathrm{i}^{(k)}</math> denote the multicomplex imaginary units; <math>\mathrm{i}^{(1)} \equiv \mathrm{i}</math>. The <math>\mathcal{C}^{(n)}_k</math> operator extracts the <math>k</math>th component of a multicomplex number of level <math>n</math>, e.g., <math>\mathcal{C}^{(n)}_0</math> extracts the real component and <math>\mathcal{C}^{(n)}_{n^2-1}</math> extracts the last, “most imaginary” component. The method can be applied to mixed derivatives, e.g. for a second-order derivative
<math display="block">\frac{\partial^2 f(x, y)}{\partial x \,\partial y} \approx \frac{\mathcal{C}^{(2)}_3(f(x + \mathrm{i}^{(1)} h, y + \mathrm{i}^{(2)} h))}{h^2}</math>

A C++ implementation of multicomplex arithmetics is available.<ref>{{cite web | last1 = Bell | first1 =  I. H.| title = mcx (multicomplex algebra library) | website = [[GitHub]]| year = 2019 | url = https://github.com/ianhbell/mcx}}</ref>

In general, derivatives of any order can be calculated using [[Cauchy's integral formula]]:<ref>Ablowitz, M. J., Fokas, A. S.,(2003). Complex variables: introduction and applications. [[Cambridge University Press]]. Check theorem 2.6.2</ref>
<math display="block">f^{(n)}(a) = \frac{n!}{2\pi i} \oint_\gamma \frac{f(z)}{(z - a)^{n+1}} \,\mathrm{d}z,</math>
where the integration is done [[Numerical integration|numerically]].

Using complex variables for numerical differentiation was started by Lyness and Moler in 1967.<ref name=LynessMoler1>{{cite journal | first1 = J. N. | last1 = Lyness | first2 = C. B. | last2 = Moler | title = Numerical differentiation of analytic functions | journal = SIAM J. Numer. Anal. | volume = 4 | year = 1967 | issue = 2 | pages = 202–210 | doi=10.1137/0704019| bibcode = 1967SJNA....4..202L }}</ref> 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.<ref>{{cite journal | title = A New Method for Generating Power Series Expansions of Functions | first1 = J | last1 = Abate | first2 = H | last2 = Dubner | journal = SIAM J. Numer. Anal. | volume =5 | issue = 1 | pages = 102–112 |date=March 1968 | doi = 10.1137/0705008| bibcode = 1968SJNA....5..102A }}</ref> An algorithm that can be used without requiring knowledge about the method or the character of the function was developed by Fornberg.<ref name=Fornberg1/>

==Differential quadrature==
Differential quadrature is the approximation of derivatives by using weighted sums of function values.<ref>Differential Quadrature and Its Application in Engineering: Engineering Applications, Chang Shu, Springer, 2000, {{isbn|978-1-85233-209-9}}.</ref><ref>Advanced Differential Quadrature Methods, Yingyan Zhang, CRC Press, 2009, {{isbn|978-1-4200-8248-7}}.</ref> 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.<ref>{{Cite journal | last1=Ahnert|first1=Karsten | last2=Abel|first2=Markus | title=Numerical differentiation of experimental data: local versus global methods | journal=Computer Physics Communications | year=2007 | volume=177 | issue=10 | pages=764–774 | doi=10.1016/j.cpc.2007.03.009 | bibcode=2007CoPhC.177..764A |s2cid=15129086 | issn=0010-4655| citeseerx=10.1.1.752.3843 }}</ref>

==See also==
*{{annotated link|Automatic differentiation}}
*{{annotated link|Five-point stencil}}
*[[List of numerical-analysis software]]
*{{annotated link|Numerical integration}}
*{{annotated link|Numerical methods for ordinary differential equations}}
*{{annotated link|Savitzky–Golay filter}}

==References==
{{reflist}}

== External links ==
{{wikibooks|Numerical Methods}}
* [http://mathworld.wolfram.com/NumericalDifferentiation.html Numerical Differentiation] from wolfram.com
* [http://www.nag.co.uk/numeric/fl/nagdoc_fl24/html/D04/d04conts.html NAG Library numerical differentiation routines]
* [http://www.boost.org/doc/libs/release/libs/math/doc/html/math_toolkit/diff.html Boost. Math numerical differentiation, including finite differencing and the complex step derivative]
*[https://sinews.siam.org/Details-Page/differentiation-without-a-difference Differentiation With(out) a Difference] by [[Nicholas Higham]], [[SIAM]] News.

{{DEFAULTSORT:Numerical Differentiation}}
[[Category:Numerical analysis]]
[[Category:Differential calculus]]