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.Example: Derivation of Central Difference
Example: ICD
\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. 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! 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.Example: Reading/Constructing the Table
pow function by just multiplying a re-used variable by 4 when calculating the extrapolations.Example: Applying Richardson’s
Example: Derivation/Reasoning