Numerical Differentiation¶
We will roughly follow the treatment in TAK § 3.2.
Goal: Given a function $f:{\mathbb R} \to {\mathbb R}$, how to you approximate its derivative?
# Standard Imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv
Loss of precision in subtracting similar numbers¶
When you subtract numbers that agree to $n$ significant digits of precision, it leads to a loss of significant digits. For example, the following two numbers have $8$ significant digits, but agree to the first three. When we subtract, the resulting number only has $5=8-3$ significant digits: $$1642.6234 - 1641.6217 = 1.0017.$$ The same of course happens for bits of precision, and for addition when adding two numbers that are nearly negatives of one another.
This issue is called loss of significance.
The definition of the derivative¶
Recall that the derivative of a function $f:{\mathbb R} \to {\mathbb R}$ at $a \in {\mathbb R}$ is $$f'(a) = \lim_{h \to 0} \frac{f(a+h)-f(a)}{h}.$$ We call $f$ differentiable at $a$ if this limit exists.
In these notes, you should think of $f$ as a mystery function (since if $f$ is given as a mathematical function, you can probably differentiate it using standard rules).
Example: Naive Computation¶
Consider the function $f(x) = e^x \sin x$ at the point $x=2.2$. We know from usual differentiation rules that $$f'(x)=e^x (\sin x+\cos x)$$ and so $f'(2.2) \approx 1.9854604310541828$ as computed below.
der = m.exp(2.2) * ( m.sin(2.2) + m.cos(2.2) )
der
But suppose we didn't know what our function was. We could attempt to compute the derivative using the approximation $$f'(0) \approx \frac{f(h)-f(0)}{h} \quad \text{for $h$ small.}$$ But how small should we make $h$? We will consider what values we get as approximate derivatives when $h$ has the form $h=2^{-k}$.
def f(x):
return np.exp(x) * np.sin(x)
Below we show a plot of $f$ near $x=0$. The plot is illustrated with the secant to the graph, a segment joining $\big(0,f(0)\big)$ to $\big(0.5, f(0.5)\big)$ whose slope is the first approximation to the derivative, $\frac{f(0.5)-f(0)}{0.5}$.
x = np.linspace(1.6, 2.8, 1001)
plt.plot(x, f(x))
plt.plot([2.2, 2.7], [f(2.2), f(2.7)])
plt.show()
Now we will form a table with our $k$ determining $h=2^{-k}$ and our approximation of $f'(0)$. Recall that $f'(0)=1$.
a = 2.2
for k in range(55):
h = 2**-k
approx = ( f(a+h) - f(a) ) / h
error = der - approx
print("k = {:>2} approx = {:>19} error = {:>10.3e}".format(k, approx, error))
Observe that the most accurate choice of $h$ is $h=2^{-27}$, and we get the derivative correct to only about $7$ or $8$ significant digits. Below, we see this error in hex, and see this translates to $25$ bits of precision, which is roughly half of the $53$ bits of precision we might expect in a float. (You can convert from number of significant digits to number of significant bits by multiplying by $\log_2 10 \approx 3.3$.)
a=2.2
k = 27
h = 2**-k
approx = ( f(a+h) - f(a) ) / h
error = der - approx
error.hex()
First we should discuss why values of $k$ far away from $27$ are so bad:
- If $k$ is small, then the secant line just isn't very close to the graph. This is illustrated in the plot above, where the slope of the secant line doesn't even have the correct sign.
- If $k$ is large, then we suffer from a loss of significance due to the subtraction in our formula.
Analyzing Loss of Significance
We will now analyze the error due to loss of significance. We will assume that near $a$, the function $f$ has the property that if $a$ and $a+h$ agree to $n$ significant bits, then $f(a)$ and $f(a+n)$ also agre to $n$ significant bits. This is the case for most well-written functions (but it depends on how a function is implemented on the computer).
Since $a=2.2$, the most significant bit has place value $2^1$. Since $h=2^{-k}$, the bits in place values $2^1, 2^0, \ldots, 2^{1-k}$ all agree between $a$ and $a+h$. Thus $a$ and $a+h$ agree to $k+1$ significant bits. Then by assumption $f(a)$ and $f(a+h)$ agree to $k+1$ significant bits. Thus computing $f(a+h)-f(a)$ results in a loss of $k+1$ significant bits. Since we are working with floats, we expect to have $53 - (k+1)= 52-k$ significant bits remaining. Division by $h$ will not result in any loss of precision (since $h$ is a power of two), so our estimate of $f'(a)$ should be correct to $52-k$ significant bits. Since $1 \leq f'(a)<2$, the place value $2^0$ is the most significant bit. So, we can only expect to be accurate to the bit representing place value $2^{k-51}$. Which means that our expected error should be the same in magnitude as $2^{k-52}$.
For large values of $k$, this explains the error. See the table below. Numbers are written in exponential notation, and observe the exponents roughly match.
a = 2.2
for k in range(30,55):
h = 2**-k
approx = ( f(a+h) - f(a) ) / h
error = der - approx
expected_error = 2**(k-52)
print("k = {:>2} error = {:>10.3e} expected error = {:.3e}". \
format(k, error, expected_error))
The list of numbers is a bit hard to understand. Below, we create a plot sending these values of $k$ to the ratio of the error and our expected error $2^{k-52}$. Observe that this ratio doesn't go much higher than $4$ in absolute value.
a = 2.2
ks = []
error_ratios = []
for k in range(30,55):
h = 2**-k
approx = ( f(a+h) - f(a) ) / h
error = der - approx
expected_error = 2**(k-52)
ks.append(k)
error_ratios.append(error/expected_error)
plt.plot(ks, error_ratios, 'o-')
plt.xlabel('$k$') # The dollar signs denote LaTeX.
plt.ylabel('Ratio of actual error to expected error')
plt.show()
So, $2^{k-52}$ is a pretty good estimate of the error. It is useful to express this in terms of $h=2^{-k}$. We have that $$\text{error} = 2^{k-52} = \frac{2^{-52}}{h}.$$
How good is a secant approximation?¶
The key to understanding how good a secant approximation to the derivative is Taylor's Theorem in the degree one case. Assuming our function is twice differentiable on the interval with endpoints $a$ and $a+h$, Taylor's Theorem says $$f(a+h) = f(a) + f'(a)h + \frac{f''(y)}{2}h^2$$ for some $y$. Manipulating this formula, we see that $$\frac{f(a+h)-f(a)}{h} = f'(a) + \frac{f''(y)}{2}h.$$ Assuming $f''(y)$ is nearly constant on the small interval with endpoints $a$ and $a+h$, we see that our error is roughly linear in $h$.
An improvement¶
The formula above actually suggests an improvement. Here is our formula above: $$\frac{f(a+h)-f(a)}{h} = f'(a) + \frac{f''(y_+)}{2}h,$$ where we have used $y_+$ instead of $y$. Applying the formula to $-h$ as well gives a $y_-$ between $a-h$ and $a$ such that $$\frac{f(a)-f(a-h)}{h} = f'(a) - \frac{f''(y_-)}{2}h.$$ Averaging these two formulas gives that $$\frac{f(a+h)-f(a-h)}{2h} = f'(a) + \left(\frac{f''(y_+)}{4}- \frac{f''(y_-)}{4}\right)h.$$ It may not be clear why this is better. But, recall that we are assuming that $h$ is so small that $f''$ is roughly constant on the interval from $a-h$ to $a+h$. Thus means that the difference $f''(y_+)-f''(y_+)$ is very close to zero.
A better understanding of the improvement. Consider applying the degree two Taylor's theorem. It says that $$f(a+h) = f(a) + f'(a)h + \frac{f''(a)}{2}h^2 + \frac{f^{(3)}(y_+)}{6}h^3$$ for some $y_+$ between $a$ and $a+h$. Similarly we see that $$f(a-h) = f(a) - f'(a)h + \frac{f''(a)}{2}h^2 - \frac{f^{(3)}(y_-)}{6}h^3$$ for some $y_-$ between $a-h$ and $a$. Subtracting the two equations yields: $$f(a+h)-f(a-h) = 2 f'(a)h + \frac{f^{(3)}(y_+)+f^{(3)}(y_-)}{6} h^3.$$ Simplifying yields that $$f'(a) = \frac{f(a+h)-f(a-h)}{2h} - \frac{f^{(3)}(y_+)+f^{(3)}(y_-)}{12} h^2.$$ Again assuming that the third derivative is roughly constant, we see that the error term is roughly proportional to $h^2$. Quadratic ($h^2$) is much better than linear ($h$) since when $h$ is small, $h^2$ is even smaller.
The two-point approximation of the derivative.
Assuming $f$ has a continuous $3$rd derivative in a neighborhood of $a$, we have $$f'(a) = \frac{f(a+h)-f(a-h)}{2h} + {\mathcal O}(h^2) \quad \text{as $h \to 0$.}$$Big O notation The formula above is an instance of Big O notation. The term ${\mathcal O}(h^2)$ can be replaced by an "error" term $E(h)$, which in this case must be $$E(h) = f'(a) - \frac{f(a+h)-f(a-h)}{2h}.$$ We say $E(h) \in {\mathcal O}(h^2)$ as $h \to 0$, which indicates that there are constants $M>0$ and a $\delta>0$ such that $$|E(h)| < M h^2 \quad \text{whenever $0<|h|<\delta$.}$$ The notation ${\mathcal O}(h^2)$ is convienient because it captures the concept that the error decays to zero at least as fast as a constant times $h^2$.
Below we illustrate a two point approximation of $f'(2.2)$, as the slope of the secant line joining $\big(2.2-0.5, f(2.2-0.5)\big)$ to $\big(2.2+0.5, f(2.2+0.5)\big)$.
x = np.linspace(1.6, 2.8, 1001)
plt.plot(x, f(x))
plt.plot(x, der*(x-2.2) + f(2.2))
plt.plot([1.7, 2.7], [f(1.7), f(2.7)])
plt.show()
Example: Testing the two-point approximation¶
We will test our approximation to the derivative using the same function $f(x) = e^x \sin x$ at the point $x=2.2$. Recall $f'(2.2) \approx 1.9854604310541828$.
We repeat the computation of the table using the new approximation:
a = 2.2
for k in range(55):
h = 2**-k
approx = ( f(a+h) - f(a-h) ) / (2*h)
error = der - approx
print("k = {:>2} approx = {:>19} error = {:>10.3e}".format(k, approx, error))
Above we see that the best approximation is when $k=18$ and the comptuation yeilds between $10$ and $11$ digits of precision. Here we look at the bits of precision in this case:
k = 18
h = 2**-k
approx = ( f(a+h) - f(a-h) ) / (2*h)
error = der - approx
error.hex()
We see that we get roughly $34$ bits of precision. We have lost only a third of the bits of precision with this formula.
Best choice of $h$¶
With the two-point approximation to $f'(a)$ there the same competing issues:
- Error due to bad approximation when choosing large $h$.
- Error loss of significance when choosing small $h$.
Essentially the best choice of $h$ is the one which makes the errors due to these two issues approximately equal. (To see this, note that for example if the greatest error is due to bad approximation, then we should shrink $h$ to make that error decrease. This will result in the error due to loss of significance to increase, but that error is smaller so it matters less.)
Even though we will consider the loss of significance due for the two-point approximation, the error estimate is the almost the same as before, roughly $2^{-52}\big/h$. From our error estimate for the two point method, we see that the error due to bad approximation is roughly $c h^2$ where $c \approx \frac{f^{(3)}(a)}{6}$. We'll just assume $c$ is a reasonably small number. So, we'd like to choose an $h$ such that roughly $$\frac{2^{-52}}{h} = c h^2.$$ Thus, $h^3 = \frac{2^{-52}}{c}$ and so $$h=\frac{2^{-17 \frac{1}{3}}}{\sqrt[3]{c}}.$$ Note that we are not interested in getting an exact answer here. We are happy simply to know what magnitude to choose for $h$. Assuming $c$ is a reasonable number, we should choose something on the order of $h=2^{-17 \frac{1}{3}}$. This will then result in roughly equal error loss to both bad approximation and loss of significance as desired.
The value $h=2^{-17 \frac{1}{3}}$ is pretty close to the $h=2^{-18}$ that we saw worked best when looking at the table.
The second derivative:¶
We will now explain how to use Taylor polynomials to estimate the 2nd derivative using three points. Evaluating the degree $2$ Taylor polynomial and considering the error gives that $$f(a+h) = f(a) + f'(a)h + \frac{f''(a)}{2}h^2 + \frac{f^{(3)}(y_+)}{6}h^3$$ for some $y_+$ between $a$ and $a+h$, and $$f(a-h) = f(a) - f'(a)h + \frac{f''(a)}{2}h^2 - \frac{f^{(3)}(y_-)}{6}h^3$$ for some $y_-$ between $a-h$ and $a$.
Adding the formulas yields: $$f(a+h) + f(a-h) = 2 f(a) + f''(a) h^2 + \frac{f^{(3)}(y_+)-f^{(3)}(y_-)}{6}h^3.$$ Solving for $f''(a)$, we see that $$f''(a) = \frac{f(a+h) - 2 f(a) + f(a-h)}{h^2} - \frac{f^{(3)}(y_+)-f^{(3)}(y_-)}{6}h.$$
This formula looks like it an error of order $h$, but if $f$ is $4$-times differentiable, we can apply the Mean value theorem to conclude that $$\frac{f^{(3)}(y_+)-f^{(3)}(y_-)}{y_+-y_-} = f^{(4)}(z)$$ for some $z$ between $y_-$ and $y_+$. Thus, $$|f^{(3)}(y_+)-f^{(3)}(y_-)| = (y_+-y_-) |f^{(4)}(z)| \leq 2h |f^{(4)}(z)|.$$ Plugging in, we see that $$\left|f''(a) - \frac{f(a+h) + 2 f(a) - f(a-h)}{h^2}\right| \leq \frac{|f^{(4)}(z)|}{3} h^2.$$ for some $z \in (-h,h)$. This is summarized below:
The three-point approximation of the 2nd derivative.
Assuming $f$ has a continuous $4$th derivative in a neighborhood of $a$, we have $$f''(a) = \frac{f(a+h) - 2 f(a) + f(a-h)}{h^2} + {\mathcal O}(h^2) \quad \text{as $h \to 0$.}$$Example
Lets return to the same example, this time to compute $f''(2.2)$. To get an exact computation, we can use calculus. Recall $$f'(x)=e^x (\sin x+\cos x).$$ Thus, $$f''(x)=e^x (\sin x+\cos x) + e^x(\cos x - \sin x)=2 e^x \cos x.$$
a = 2.2
der2 = 2 * m.exp(a) * m.cos(a)
der2
Here is a table computing approximates for $f''(a)$ using for values of $h=2^{-k}$. Note that the smallest error we see is about $10^{-7}$ when $k=13$ or so.
a = 2.2
for k in range(26):
h = 2**-k
approx = (f(a+h) + f(a-h) - 2*f(a)) / h**2
error = der2 - approx
print("k = {:>2} approx = {:>19} error = {:>10.3e}".format(k, approx, error))
We can get a better approximation using more points as indicated below.
The five-point approximation of the 2nd derivative.
Assuming $f$ has a continuous $6$th derivative in a neighborhood of $a$, we have $$f''(a) = \frac{-f(a+2h) + 16 f(a+h) - 30 f(a) + 16 f(a-h) - f(a-2h)}{12 h^2} + {\mathcal O}(h^4) \quad \text{as $h \to 0$.}$$Proof:
I will briefly explain how this is obtained. Considering Taylor's theorem for $f(a \pm h)$ we can see that $$f(a+h)+f(a-h)-2 f(a) = f''(a)h^2 + \frac{f^{(4)}(a)}{12} h^4 + \frac{f^{(5)}(y_1) - f^{(5)}(y_{-1})}{5!}h^5,$$ for some $y_1 \in (a,a+h)$ and $y_{-1} \in (a-h,a)$. Applying the Mean value theorem then tell us that $$f^{(5)}(y_1) - f^{(5)}(y_{-1})=(y_1 - y_{-1})f^{(6)}(z_1)$$ for some $z_1 \in (y_{-1}, y_1).$ Then, $$f(a+h)+f(a-h)-2 f(a) = f''(a)h^2 + \frac{f^{(4)}(a)}{12} h^4 + \frac{(y_1 - y_{-1})f^{(6)}(z_1)}{5!}h^5.$$
Then substituting $2h$ for $h$, we see that $$f(a+2h)+f(a-2h)-2 f(a) = 4 f''(a)h^2 + \frac{4 f^{(4)}(a)}{3} h^4 + \frac{32(y_2 - y_{-2})f^{(6)}(z_2)}{5!}h^5,$$ for some $y_2 \in (a,a+2h)$, $y_{-2} \in (a-2h,a)$, and $z_2 \in (y_{-2},y_2)$.
Letting $X_1=f(a+h)+f(a-h)-2 f(a)$ and $X_2 = f(a+2h)+f(a-2h)-2 f(a)$, we see that we can remove the degree four term as follows:
$$16 X_1 - X_2 = 12 f''(a)h^2 + \frac{16(y_1 - y_{-1})f^{(6)}(z_1)-32(y_2 - y_{-2})f^{(6)}(z_2)}{5!}h^5.$$
Solving for $f''(a)$, we see that
$$f''(a) = \frac{16 X_1 - X_2}{12 h^2} - \frac{(y_1 - y_{-1})f^{(6)}(z_1)-2(y_2 - y_{-2})f^{(6)}(z_2)}{90}h^3.$$
Now observe that $0
Example
Here is a table computing approximates for $f''(a)$ using the five point formula for values of $h=2^{-k}$. Note that the smallest error we see is about $10^{-10}$ when $k=9$ or so. This is a bit better than the three point fomula.
a = 2.2
for k in range(26):
h = 2**-k
approx = (-f(a+2*h) + 16*f(a+h) -30 *f(a) + 16*f(a-h) - f(a-2*h)) / (12*h**2)
error = der2 - approx
print("k = {:>2} approx = {:>19} error = {:>10.3e}".format(k, approx, error))