Goal: Compute the numerical derivative of a function f(x) at a particular point (x_0).
1. Initialize
Pick an initial step size h and tolerance \epsilon
2. Evaluate estimation formula
\boxed{ f'(x_0)_i = \frac{f(x_0 + h/2) - f(x_0 - h/2)}{h} } \\~\\ \small\textit{where $i$ current num. of iterations}
3. Stopping criteria
If \lvert f'_i(x_0) - f'_{i-1} (x_0) \rvert < \epsilon, stop.
Otherwise, h \leftarrow h/2 and go to step 2.
We’re going to start with an approximation of f’(x) we previously saw used by the Secant method.
f'(x) \approx \frac{f(x+h) - f(x)}{h}
This formula is one-sided (it only approaches from the right). To improve accuracy, we can sandwich the estimation by splitting h in half, to approachg from both sides.
The resulting sandwiching formula is:
f'(x_0) = \frac{f(x_0 + h/2) - f(x_0 - h/2)}{h}
Note — This also ends just being more accurate than one-sided because it cancels out more error terms.
The smaller h is, the better the approximation. Our stopping criteria will be when the delta between iterations falls below a tolerance \epsilon.
Q:
f(x) = x^2 + \sin x, x_0 = 1.7, \epsilon = 0.00001
Note — Professor picked this polynomial because smooth polynomials converge really fast.
Analytic Derivative — Just so we have something to compare against, f'(1.7) = 3.2711555....
A:
It 1: h = 0.5
\begin{aligned} f'_1 (x_0) &= \frac{ f(1.7 + 0.25) - f(1.7 - 0.25) }{ 0.5 } \\ &= 3.27493448... \end{aligned}
It 2: h = 0.25
\begin{aligned} f'_2 (x_0) &= \frac{ f(1.7 + 0.125) - f(1.7 - 0.125) }{ 0.25 } \\ &= 3.271490776... \end{aligned}
\begin{aligned} \Delta &= |3.271490776 - 3.27493448| \\ &= 0.001002672 \not\lt \epsilon \end{aligned}
It 3: h = 0.125
\begin{aligned} f'_3 (x_0) &= \frac{ f(1.7 + 0.0625) - f(1.7 - 0.0625) }{ 0.125 } \\ &= 3.271239372... \end{aligned}
\begin{aligned} \Delta &= |3.271239372 - 3.271490776| \\ &= 0.000251... \not\lt \epsilon \end{aligned}
It 4: h = 0.0625
\begin{aligned} f'_4 (x_0) &= \frac{ f(1.7 + 0.03125) - f(1.7 - 0.03125) }{ 0.0625 } \\ &= 3.271176475... \end{aligned}
\begin{aligned} \Delta &= |3.271176475 - 3.271239372| \\ &= 0.000062897... \not\lt \epsilon \end{aligned}
It 5: h = 0.03125
\begin{aligned} f'_5 (x_0) &= \frac{ f(1.7 + 0.015625) - f(1.7 - 0.015625) }{ 0.03125 } \\ &= 3.27116074834... \end{aligned}
\begin{aligned} \Delta &= |3.27116074834 - 3.271176475| \\ &= 0.00001572... \not\lt \epsilon \end{aligned}
We’ll stop here, but this should get the idea across. Pretty simple, but generally slow to converge.
\boxed{ \text{Richardson's} \begin{cases} D(i,0) &= \frac{ f(x_0 + k) - f(x_0 - k)}{2k} \\ D(i,j) &= D(i,j-1) + \frac{D(i, j-1) - D(i-1, j-1)}{4^j - 1} \end{cases} } \\ \small\textit{where initial step size: $h$; step size for row $i$: $k = h/2^i$}
Algorithm:
\begin{array}{l | c c c c l} & \mathbf{j=0} & \mathbf{j=1} & \mathbf{j=2} & \cdots & \mathbf{j=N} \\ \hline \mathbf{i=0} & D(0,0) \\ \mathbf{i=1} & D(1,0) & D(1,1) \\ \mathbf{i=2} & D(2,0) & D(2,1) & D(2,2) \\ \mathbf{i=3} & D(3,0) & D(3,1) & D(3,2) & \ddots \\ \vdots & \vdots & \vdots & \ddots & \ddots \\ \mathbf{i=N} & D(N,0) & D(N,1) & \cdots & \cdots & D(N,N) \end{array}
Column zero (j=0) is the base calculation, which is just the previous central difference formula.
The other columns (j>0) are the extrapolations. Each new D(i,j) depends on the value to its immediate left (D(i,j-1)) and diagonally up-left (D(i-1,j-1)).
The diagonal elements (D(i,i)) are the most-accurate approximations at each iteration i.
Aside — Overall, the table is just to make calculations on-paper easier. Have fun writing a bajillion numbers regardless.
- When programming, note how you (1) only need to store the two most-recent rows of the table and (2) how the 4^j term can be calculated without repeatedly calling a
powfunction by just multiplying a re-used variable by 4 when calculating the extrapolations.
Q:
f(x) = x^2 + \sin x, x_0 = 1.7, \epsilon = 0.00001
h = 0.25
Analytic Derivative — Just so we have something to compare against, f'(1.7) = 3.2711555....
A:
Round 1
\begin{aligned} D(0,0) &= \frac{ f(x_0 + h_0) - f(x_0 - h_0) }{ 2h_0 } \\ &= \frac{ f(1.7 + 0.25) - f(1.7 - 0.25) }{ 0.5 } \\ &= 3.27493448... \end{aligned}
\begin{aligned} D(1,0) &= \frac{ f(1.7 + 0.125) - f(1.7 - 0.125) }{ 0.25 } \\ &= 3.271490776... \end{aligned}
Now we’ll do something new, the extrapolation term.
\begin{aligned} D(1,1) &= D(1,0) + \frac{ D(1,0) - D(0,0) }{ 3 } \\ &= 3.271490776 + \frac{ 3.271490776 - 3.27493448 }{ 3 } \\ &= 3.271156552... \end{aligned}
At the first round, we’ve already gotten more digits than the precious method at 5 iterations.
Round 2
D(2,0)=3.271239372
I think this was pulled form the second method?
\begin{aligned} D(2,1) &= D(2,0) + \frac{D(2,0) - D(1,0)}{3} \\ &= 3.271239372 + \frac{3.271239372 - 3.271490776}{3} \\ &= 3.271155571... \end{aligned}
Round 3:
\begin{aligned} D(2,2) &= D(2,1) + \frac{D(2,1) - D(1,1)}{15} \\ &= 3.271155571 + \frac{3.271155571 - 3.271156552}{15} \\ &= 3.271155506 \end{aligned}
It has fully converged to the analytic derivative!
- Now, we know we converged because we have the analytic derivative, but in the real-world, we wouldn’t know when to stop. For that we’d do like at most five, or go until we start getting the same number with minor fluctuations (then we just pick one).
As we iterated in our previous method, the error followed a predictable pattern.
We can take the two “wrong” answers (D(i-1,j-1) and D(i,j-1)) and combine them to get a more correct answer (D(i,j)) by cancelling-out the largest part of the error.
Skipping through the derivation, we end up with this “accelerator”:
D(i,j) = D(i,j-1) + \frac{ D(i,j-1) - D(i-1,j-1) }{ 4^j - 1 }
Each new column j in the table represents a higher level of extrapolation and converges must faster than the column before it, and the diagonal D(i,j) represents the best possible estimate at each iteration i.