4.2. Approximation of the second derivative of a function#
Approximation of the second derivative can be obtained using the forward, back, and central approximation of the first derivatives. Note that,
Second derivative central difference:
\[\begin{align*}
f''(x)&\approx {\frac {\delta _{h}^{2}[f](x)}{h^{2}}}=\frac {1}{h}\left({\frac {f(x+2h)-f(x+h)}{h}}-{\frac {f(x+h)-f(x)}{h}}\right)
\\ &={\frac {f(x+2h)-2f(x+h)+f(x)}{h^{2}}}.
\end{align*}\]
Second derivative forward difference:
\[\begin{align*}
f''(x) &\approx {\frac {\Delta _{h}^{2}[f](x)}{h^{2}}}=\frac {1}{h}\left({\frac {f(x)-f(x-h)}{h}}-{\frac {f(x-h)-f(x-2h)}{h}}\right)
\\&={\frac {f(x)-2f(x-h)+f(x-2h)}{h^{2}}}.
\end{align*}\]
Second derivative backward difference:
\[\begin{align*}
f''(x)&\approx {\frac {\nabla _{h}^{2}[f](x)}{h^{2}}}=\frac {1}{h}\left({\frac {f(x+h)-f(x)}{h}}-{\frac {f(x)-f(x-h)}{h}}\right)
\\&={\frac {f(x+h)-2f(x)+f(x-h)}{h^{2}}}.
\end{align*}\]
Example:
Consider \(f(x) = 3x\exp(x) - \cos(x) + \sin(x)\). This function is defined and continuous on [0,1]. Discretize \([0,1]\) using \(h=0.01\) and approximate \(f''(x)\) using
Forward Difference approximation of the 1st derivative at \(\{0.02, 0.03, 0.04, \ldots,0.97,0.98 \}\).
Backward Difference approximation of the 1st derivative at \(\{0.02, 0.03, 0.04, \ldots,0.97,0.98 \}\).
Central Difference approximation of the 1st derivative at \(\{0.02, 0.03, 0.04, \ldots,0.97,0.98 \}\).
Solution:
Second derivative forward difference:
\[\begin{align*}
f''(0.02) & \approx \frac {f(0.04)-2f(0.03)+f(0.02)}{(0.01)^{2}} = 7.275114,\\
f''(0.03) & \approx \frac {f(0.05)-2f(0.04)+f(0.03)}{(0.01)^{2}} = 7.369059,\\
f''(0.04) & \approx \frac {f(0.06)-2f(0.05)+f(0.04)}{(0.01)^{2}} = 7.464166,\\
&\vdots\\
f''(0.96) & \approx \frac {f(0.98)-2f(0.97)+f(0.96)}{(0.01)^{2}} = 24.069708,\\
f''(0.97) & \approx \frac {f(0.99)-2f(0.98)+f(0.97)}{(0.01)^{2}} = 24.377588 ,\\
f''(0.98) & \approx \frac {f(1.00)-2f(0.99)+f(0.98)}{(0.01)^{2}} = 24.689394.
\end{align*}\]
Similarly,
import numpy as np
import pandas as pd
# f(x)
f = lambda x: 3*x*np.exp(x) - np.cos(x)
# f''(x)
f2 = lambda x: np.cos(x) + 6*np.exp(x) + 3*x*np.exp(x)
a = 0
b = 1
h = 1e-2
xn = np.arange(a, b + h, h)
yn = f(xn)
Consider the inner points of this example. For these points
from IPython.display import display, Latex
Ind = np.arange(2, len(xn)-2)
# Exact
Exact = f2(xn[Ind])
# Forward difference
Forward = (yn[Ind+2]-2*yn[Ind+1]+yn[Ind])/(h**2)
# Backward difference
Backward = (yn[Ind]-2*yn[Ind-1]+yn[Ind-2])/(h**2)
# Central difference
Central = (yn[Ind+1]-2*yn[Ind]+yn[Ind-1])/(h**2)
display(Latex('''\\max_{2\leq i\leq N-2}\\left|f''(x_i)-\dfrac {\\Delta^2_{h}[f](x_i)}{h^2} \\right| = %.4e''' % max(abs(Forward - Exact))))
display(Latex('''\\max_{2\leq i\leq N-2}\\left|f''(x_i)-\dfrac {\\nabla^2_{h}[f](x_i)}{h^2} \\right| = %.4e''' % max(abs(Backward - Exact))))
display(Latex('''\\max_{2\leq i\leq N-2}\\left|f''(x_i)-\dfrac {\\delta^2_{h}[f](x_i)}{h^2} \\right| = %.4e''' % max(abs(Central - Exact))))
\[\max_{2\leq i\leq N-2}\left|f''(x_i)-\dfrac {\Delta^2_{h}[f](x_i)}{h^2} \right| = 3.1213e-01\]
\[\max_{2\leq i\leq N-2}\left|f''(x_i)-\dfrac {\nabla^2_{h}[f](x_i)}{h^2} \right| = 3.0755e-01\]
\[\max_{2\leq i\leq N-2}\left|f''(x_i)-\dfrac {\delta^2_{h}[f](x_i)}{h^2} \right| = 3.2708e-04\]
4.2.1. Order of Convergence#
Let \(f_{i+j}\) be \(f(x_i+ jh)\) for \(j \in \mathbb{Z}\). Using Taylor theorem, we have,
\[\begin{align*}
f_{i+2} & = f_{i} + 2hf_{i}' + h^2 f_{i}'' + \frac{2}{3} h^3 f_{i}^{(3)} + \frac{1}{3} h^4 f_{i}^{(4)} + O(h^5),\\
f_{i+1} & = f_{i} + hf_{i}' + \frac{1}{2} h^2 f_{i}'' + \frac{1}{6} h^3 f_{i}^{(3)} + \frac{1}{24} h^4 f_{i}^{(4)} + O(h^5),\\
f_{i-1} & = f_{i} - hf_{i}' + \frac{1}{2} h^2 f_{i}'' - \frac{1}{6} h^3 f_{i}^{(3)} + \frac{1}{24} h^4 f_{i}^{(4)} + O(h^5),\\
f_{i-2} & = f_{i} - 2hf_{i}' + h^2 f_{i}'' - \frac{2}{3} h^3 f_{i}^{(3)} + \frac{1}{3} h^4 f_{i}^{(4)} + O(h^5).
\end{align*}\]
and
Forward difference:
(4.22)#\[\begin{align}
f''(x_i)\approx {\frac {\Delta _{h}^{2}[f](x)}{h^{2}}}&=\frac {f_{i+2}-2f_{i+1}+f_{i}}{h^{2}} \notag
%\\ & = \frac {1}{h^{2}}\left[f_{i+2}-2f_{i+1}+f_{i}\right] \notag
= \frac {1}{h^{2}}\left[h^2 f_{i}'' + h^3 f_{i}^{(3)} + \frac{7}{2} h^4 f_{i}^{(4)} + \mathcal{O}(h^5)\right] \notag
\\ & = f_{i}'' + h f_{i}^{(3)} + \frac{7}{12} h^2 f_{i}^{(4)} + O(h^3),
\end{align}\]
Thus,
(4.23)#\[\begin{equation}
\left(f''(x_i)-\dfrac {\Delta^2_{h}[f](x_i)}{h^2} \right)=
{{h} f_{i}^{(3)}} + \frac{7}{12} h^2 f_{i}^{(4)} + \mathcal{O}(h^3),
\end{equation}\]
Backward difference
(4.24)#\[\begin{align}
f''(x)\approx {\frac {\nabla _{h}^{2}[f](x)}{h^{2}}} & =\frac {f_{i}-2f_{i-1}+f_{i-2}}{h^{2}}
= \frac {1}{h^{2}}\left[ f_{i}'' -h^3 f_{i}^{(3)} + \frac{7}{12} h^4 f_{i}^{(4)} + \mathcal{O}(h^5)\right]
\notag
\\ &
= f_{i}'' -h f_{i}^{(3)} + \frac{7}{12} h^2 f_{i}^{(4)} + \mathcal{O}(h^3),
\end{align}\]
Therefore,
(4.25)#\[\begin{equation}
\left(f''(x_i)-\dfrac {\nabla^2_{h}[f](x_i)}{h^2} \right)=
{-{h} f_{i}^{(3)}} + \frac{7}{12} h^2 f_{i}^{(4)} + \mathcal{O}(h^3),
\end{equation}\]
and,
(4.26)#\[\begin{equation}
\left|f''(x_i)-\dfrac {\nabla^2_{h}[f](x_i)}{h^2} \right| = \mathcal{O}(h).
\end{equation}\]
Central difference
(4.27)#\[\begin{align}
f''(x)\approx {\frac {\delta _{h}^{2}[f](x)}{h^{2}}}& =\frac {f_{i+1}-2f_{i}+f_{i-1}}{h^{2}} = \frac {1}{h^{2}}\left[ f_{i+1}-2f_{i}+f_{i-1}\right]
\notag
\\ & = \frac {1}{h^{2}}\left[ f_{i}'' + \frac{1}{12} h^4 f_{i}^{(4)} + \mathcal{O}(h^5)\right]
= f_{i}'' \frac{1}{12} h^2 f_{i}^{(4)} + O(h^3),
\end{align}\]
Then,
(4.28)#\[\begin{equation}
\left(f''(x_i)-\dfrac {\delta^2_{h}[f](x_i)}{h^2} \right)=
{ \frac{1}{12} {h^2} f_{i}^{(4)}} + \mathcal{O}(h^3),
\end{equation}\]
and,
(4.29)#\[\begin{equation}
\left|f''(x_i)-\dfrac {\delta^2_{h}[f](x_i)}{h^2} \right|= \mathcal{O}(h^{2}).
\end{equation}\]
To demonstrate the order of convergence numerically, we define
\[\begin{align*}
E_{h}^{F} &= \max_{i} \left|f''(x_i)-\dfrac {\Delta^2_{h}[f](x_i)}{h^2} \right|,\\
E_{h}^{B} &= \max_{i} \left|f''(x_i)-\dfrac {\nabla^2_{h}[f](x_i)}{h^2} \right|,\\
E_{h}^{C} &= \max_{i} \left|f''(x_i)-\dfrac {\delta^2_{h}[f](x_i)}{h^2} \right|.
\end{align*}\]
Then, It follows that,
(4.30)#\[\begin{align}
\begin{cases}
\dfrac{E_{h}^{F}}{E_{h/2}^{F}} \approx 2^1,\\
\dfrac{E_{h}^{B}}{E_{h/2}^{B}} \approx 2^1,\\
\dfrac{E_{h}^{C}}{E_{h/2}^{C}} \approx 2^2
\end{cases}
\quad \Rightarrow \quad
\begin{cases}
\alpha_{\text{F}} = \log_{2}\left(\dfrac{E_{h}^{F}}{E_{h/2}^{F}} \right),\\
\alpha_{\text{B}} = \log_{2}\left(\dfrac{E_{h}^{B}}{E_{h/2}^{B}} \right),\\
\alpha_{\text{C}} = \log_{2}\left(\dfrac{E_{h}^{C}}{E_{h/2}^{C}} \right).
\end{cases}
\end{align}\]
From the theoretical part, we expect to see \(\alpha_{\text{F}} \approx 1\), \(\alpha_{\text{B}} \approx 1\), and \(\alpha_{\text{C}} =2 \).
Example: For the above example, test the order of convergence of the forward difference, the backward difference, and the central difference numerically.
Solution: We create a table in the form of the following table for this Example. Because of the nature of this question, it would be simpler to try this in Python, MATLAB, or another programming language.
import numpy as np
import pandas as pd
# f(x)
f = lambda x: 3*x*np.exp(x) - np.cos(x)
# f''(x)
f2 = lambda x: np.cos(x) + 6*np.exp(x) + 3*x*np.exp(x)
a = 0; b = 1; h = 0.1; t = 0.5;
h = [2**(-i) for i in range(4, 12)]
Cols = ['h', 'EhF', 'EhB', 'EhC']
Table = pd.DataFrame(np.zeros([len(h), len(Cols)], dtype = float), columns=Cols)
Table[Cols[0]] = h
for i in range(len(h)):
xn = np.arange(a, b, h[i])
yn = f(xn)
Ind = np.arange(2, len(xn)-2)
# Exact
Exact = f2(xn[Ind])
# Forward difference
Forward = (yn[Ind+2]-2*yn[Ind+1]+yn[Ind])/(h[i]**2)
Table[Cols[1]][i] = max(abs(Forward - Exact))
# Backward difference
Backward = (yn[Ind]-2*yn[Ind-1]+yn[Ind-2])/(h[i]**2)
Table[Cols[2]][i] = max(abs(Backward - Exact))
# Central difference
Central = (yn[Ind+1]-2*yn[Ind]+yn[Ind-1])/(h[i]**2)
Table[Cols[3]][i] = max(abs(Central - Exact))
del Forward, Backward, Central, Exact, xn, yn, Ind
Table.insert(1, 'N', ((b-a)/Table['h']).astype(int))
display(Table.style.set_properties(subset=['h', 'N'], **{'background-color': 'PaleGreen', 'color': 'Black',
'border-color': 'DarkGreen'}).format(dict(zip(Table.columns.tolist()[-3:], 3*["{:.4e}"]))))
| h | N | EhF | EhB | EhC | |
|---|---|---|---|---|---|
| 0 | 0.062500 | 16 | 1.6406e+00 | 1.4954e+00 | 1.0369e-02 |
| 1 | 0.031250 | 32 | 9.0253e-01 | 8.6172e-01 | 2.9146e-03 |
| 2 | 0.015625 | 64 | 4.7335e-01 | 4.6254e-01 | 7.7237e-04 |
| 3 | 0.007812 | 128 | 2.4240e-01 | 2.3962e-01 | 1.9879e-04 |
| 4 | 0.003906 | 256 | 1.2266e-01 | 1.2195e-01 | 5.0424e-05 |
| 5 | 0.001953 | 512 | 6.1696e-02 | 6.1519e-02 | 1.2698e-05 |
| 6 | 0.000977 | 1024 | 3.0941e-02 | 3.0896e-02 | 3.1860e-06 |
| 7 | 0.000488 | 2048 | 1.5493e-02 | 1.5482e-02 | 8.0559e-07 |
# This part is used for producing tables and figures
import sys
sys.path.insert(0,'..')
import hd_tools as hd
from bokeh.plotting import show
hd.derivative_ConvergenceOrder(vecs = [Table['EhF'].values, Table['EhB'].values, Table['EhC'].values],
labels = ['Forward Difference', 'Backward Difference', 'Central Difference'],
xlabel = r"$$i$$",
ylabel = r"$$\log_{2} \left( E_{h_{i}} / E_{h_{i-1}} \right)$$",
title = 'Order of accuracy: Forward, Backward and Central Differences',
legend_orientation = 'horizontal', ylim = [0, 2.1])