Editing Numerical differentiation (section)

==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>

=== Numerical Example ===
Consider approximating the derivative of <math>f(x)=x \sin{x}</math> at the point <math>x_{0} = \frac{\pi}{4}</math>.
Since <math>f'(x)=\sin{x} + x \cos{x}</math>, the exact value is
<math display="block">
f'(\frac{\pi}{4}) = \sin{\frac{\pi}{4}} + \frac{\pi}{4}\cos{\frac{\pi}{4}} = \frac{1}{\sqrt{2}} + \frac{\pi}{4\sqrt{2}} \approx 1.2624671484563432.
</math>

{| class=wikitable style="border: none;"
! scope=col | Formula
! scope=col | h
! scope=col | Approximation
! scope=col | Error
|-
| rowspan=4 | Three-point forward difference formula || 0.1 || 1.2719084899816118 || 9.441 x 10<sup>-3</sup>
|-
| 0.01 || 1.2625569346253918 || 8.978 x 10<sup>-5</sup>
|-
| 0.001 || 1.2624680412510747 || 8.927 x 10<sup>-7</sup>
|-
| 0.0001 || 1.2624671573796542 || 8.923 x 10<sup>-9</sup>
|-
| rowspan=4 | Three-point backward difference formula || 0.1 || 1.2580094219247624 || 4.457 x 10<sup>-3</sup>
|-
| 0.01 || 1.2624225374520737 || 4.461 x 10<sup>-5</sup>
|-
| 0.001 || 1.2624667023429792 || 4.461 x 10<sup>-7</sup>
|-
| 0.0001 || 1.2624671439953605 || 4.460 x 10<sup>-9</sup>
|-
| rowspan=4 | Three-point central difference formula || 0.1 || 1.2707750261498707 || 8.307 x 10<sup>-3</sup>
|-
| 0.01 || 1.2625557981227442 || 8.864 x 10<sup>-5</sup>
|-
| 0.001 || 1.2624680401146504 || 8.916 x 10<sup>-7</sup>
|-
| 0.0001 || 1.262467157379099 || 8.922 x 10<sup>-9</sup>
|}

=== Code ===
The following is an example of a [[Python (programming language)|Python]] implementation for finding derivatives numerically for <math>f(x) = \frac{2x}{1+\sqrt{x}}</math> using the various three-point difference formulas at <math>x_{0} = 4</math>. The function <code>func</code> has derivative <code>func_prime</code>.

{| role="presentation" class="wikitable mw-collapsible mw-collapsed"
|'''Example implementation in [[Python (programming language)|Python]]'''
|-
|<syntaxhighlight lang="python3" line="1">
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_backward(value, h):
    return ((-1/2)*func(value - h) + (1/2)*func(value + h)) / h

def three_point_central(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("============================================")
</syntaxhighlight>
|}

{| role="presentation" class="wikitable mw-collapsible mw-collapsed"
|'''Output'''
|-
|<syntaxhighlight lang="python3" line="1">
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
============================================
</syntaxhighlight>
|}

=== 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>
-->