4.3. Higher-order approximations#
We’ve previously covered forward, backward, and central difference approximations for a function’s first and second derivatives. These approximations have convergence orders of 1 and 2, but what about approximations with greater convergence rates?
In this section, we discuss a general method for developing these methods. Recall Taylor Polynomials from Chapter 1.
Assume \(f \in C^n[a, b]\) and let \(f^{(n)}\) denotes \(n\)th derivative of the function \(f\). If \(f^{(n)}\) also exists on \([a,b]\), then for any \(x\) from the interval \([a,b]\), there exists \(c \leq \xi(x) \leq x\) such that
(4.31)#\[\begin{align}
f(x) = \underbrace{\sum_{j = 0}^{n} \frac{f^{(j)}(c)}{j!} (x - c)^{j} }_{P_{n}(x):~\text{nth Taylor polynomial}} +
\underbrace{\frac{f^{(n+1)} (\xi(x))}{(n+1)!} {(x - c)^{n+1}}}_{R_{n}(x):~\text{truncation error}}
\end{align}\]
Taking the limit of \(P_{n}(x)\) as \(n \to \infty\) gives an infinite series which is known as \textbf{Taylor series} for \(f\) about \(c\).
These approximations are with convergence order of 1 and 2; however what about higher convergne rate approximations.
4.3.1. Higher-order approximations of the first derivative of a function#
Now, it follows from \(x = x_{i} + h\) and \(c = x_{i}\) that,
(4.32)#\[\begin{align}
f(x_{i} + h) & = \sum_{j = 0}^{n} \frac{f^{(j)}(x_{i})}{j!} (x_{i} + h - x_{i})^{j}
+ \mathcal{O}(h^{n+1})
\notag \\ &
= \sum_{j = 0}^{n} \frac{f^{(j)}(x_{i})}{j!} h^{j} + \mathcal{O}(h^{n+1}).
\end{align}\]
Similarly,
(4.33)#\[\begin{align}
f(x_{i} - 2\,h) & = \sum_{j = 0}^{n} \frac{f^{(j)}(x_{i})}{j!} (-2\,h)^{j} + \mathcal{O}(h^{n+1}),\\
f(x_{i} - h) & = \sum_{j = 0}^{n} \frac{f^{(j)}(x_{i})}{j!} (-h)^{j} + \mathcal{O}(h^{n+1}),\\
f(x_{i} + h) & = \sum_{j = 0}^{n} \frac{f^{(j)}(x_{i})}{j!} (h)^{j} + \mathcal{O}(h^{n+1}),\\
f(x_{i} + 2\,h) & = \sum_{j = 0}^{n} \frac{f^{(j)}(x_{i})}{j!} (2\,h)^{j} + \mathcal{O}(h^{n+1}).
\end{align}\]
Let \(f_{i+j}\) be \(f(x_i+ jh)\) for \(j \in \mathbb{Z}\). The above can be polynomials can be expressed as follows,
(4.34)#\[\begin{align}
\begin{cases}
f_{i+2} & = f_{i} + 2hf_{i}' + 2 h^2 f_{i}'' + \dfrac{4}{3} h^3 f_{i}^{(3)} + \dfrac{2}{3} h^4 f_{i}^{(4)} + \mathcal{O}(h^5),
\\
f_{i+1} & = f_{i} + hf_{i}' + \dfrac{1}{2} h^2 f_{i}'' + \dfrac{1}{6} h^3 f_{i}^{(3)} + \dfrac{1}{24} h^4 f_{i}^{(4)} + \mathcal{O}(h^5),
\\
f_{i-1} & = f_{i} - hf_{i}' + \dfrac{1}{2} h^2 f_{i}'' - \dfrac{1}{6} h^3 f_{i}^{(3)} + \dfrac{1}{24} h^4 f_{i}^{(4)} + \mathcal{O}(h^5),
\\
f_{i-2} & = f_{i} - 2hf_{i}' + 2 h^2 f_{i}'' - \dfrac{4}{3} h^3 f_{i}^{(3)} + \dfrac{2}{3} h^4 f_{i}^{(4)} + \mathcal{O}(h^5). \end{cases}
\end{align}\]
For example, to develop a fourth-order approximation for the first derivative of a function using five points, we can solve the following linear system for \(a_{j}\) with \(0\leq j \leq 4\).
(4.35)#\[\begin{align}\label{4th-order-eq}
h f'(x_i) = \sum_{j = -2}^{2} a_{j+2} f_{i+j} & = a_{0}f_{i-2} + a_{1}f_{i-1} + a_{2}f_{i} + a_{3}f_{i+1} + a_{4}f_{i+2}
\notag \\ &
= a_{0}\left(f_{i} - 2hf_{i}' + 2h^2 f_{i}'' - \frac{4}{3} h^3 f_{i}^{(3)} + \frac{2}{3} h^4 f_{i}^{(4)} \right)
\notag \\ &
+ a_{1}\left( 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)} \right)
\notag \\ &
+ a_{2}f_{i}
\notag \\ &
+ a_{3}\left( 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)} \right)
\notag \\ &
+ a_{4}\left( f_{i} + 2hf_{i}' + 2h^2 f_{i}'' + \frac{4}{3} h^3 f_{i}^{(3)} + \frac{2}{3} h^4 f_{i}^{(4)} \right)+ \mathcal{O}(h^5).
\end{align}\]
In other words,
(4.36)#\[\begin{align}
h f'(x_i) &= \left(a_{0}f_{i} - 2h\,a_{0}f_{i}' + 2h^2\,a_{0} f_{i}'' - \frac{4}{3} h^3 \,a_{0}f_{i}^{(3)} + \frac{2}{3} h^4 \,a_{0}f_{i}^{(4)} \right)
\notag \\ &
+ \left( a_{1}f_{i} - h\,a_{1}f_{i}' + \frac{1}{2} h^2 \,a_{1} f_{i}'' - \frac{1}{6} h^3 \,a_{1} f_{i}^{(3)} + \frac{1}{24} h^4 \,a_{1} f_{i}^{(4)} \right)
\notag \\ &
+ a_{2}f_{i}
\notag \\ &
+ \left( a_{3} f_{i} + h\,a_{3}f_{i}' + \frac{1}{2} h^2 \,a_{3} f_{i}'' + \frac{1}{6} h^3 \,a_{3} f_{i}^{(3)} + \frac{1}{24} h^4 \,a_{3} f_{i}^{(4)} \right)
\notag \\ &
+ \left( a_{4}f_{i} + 2h\,a_{4}f_{i}' + 2h^2 \,a_{4}f_{i}'' + \frac{4}{3} h^3 \,a_{4}f_{i}^{(3)} + \frac{2}{3} h^4 \,a_{4}f_{i}^{(4)} \right)+ \mathcal{O}(h^5).
\end{align}\]
or
(4.37)#\[\begin{align}
h f'(x_i) &= \left(a_{0} + a_{1} + a_{2} + a_{3} + a_{4}\right)f_{i}
\notag \\ &
+ \left(- 2a_{0} - a_{1} + 0\,a_{2} + a_{3} + 2a_{4} \right)h\,f_{i}'
\notag \\ &
+ \left(2a_{0} +\dfrac{1}{2} a_{1} + 0a_{2} + \dfrac{1}{2} a_{3} + 2a_{4}\right)h^2 \,f_{i}''
\notag \\ &
+ \left(-\dfrac{4}{3}a_{0} -\dfrac{1}{6} a_{1} + 0a_{2} + \dfrac{1}{6} a_{3} + \dfrac{4}{3}a_{4} \right)h^3 \,f_{i}^{(3)}
\notag \\ &
+\left(\dfrac{2}{3}a_{0} +\dfrac{1}{24} a_{1} + 0a_{2} + \dfrac{1}{24} a_{3} + \dfrac{2}{3}a_{4}\right)h^4 \,f_{i}^{(4)} + \mathcal{O}(h^5).
\end{align}\]
Now since \(h>0\), it follows that
(4.38)#\[\begin{align}
\begin{cases}
a_{0} + a_{1} + a_{2} + a_{3} +a_{4} &=0, \\
-2a_{0} - a_{1} + 0a_{2} + a_{3} +2a_{4} &=1, \\
2a_{0} +\dfrac{1}{2} a_{1} + 0a_{2} + \dfrac{1}{2} a_{3} + 2a_{4} &=0, \\
-\dfrac{4}{3}a_{0} -\dfrac{1}{6} a_{1} + 0a_{2} + \dfrac{1}{6} a_{3} + \dfrac{4}{3}a_{4} &=0, \\
\dfrac{2}{3}a_{0} +\dfrac{1}{24} a_{1} + 0a_{2} + \dfrac{1}{24} a_{3} + \dfrac{2}{3}a_{4} &=0.
\end{cases}
\end{align}\]
The solution to the system can be found as follows,
(4.39)#\[\begin{equation}
\left[\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1\\ -2 & -1 & 0 & 1 & 2\\ 2 & \frac{1}{2} & 0 & \frac{1}{2} & 2\\ -\frac{4}{3} & -\frac{1}{6} & 0 & \frac{1}{6} & \frac{4}{3}\\ \frac{2}{3} & \frac{1}{24} & 0 & \frac{1}{24} & \frac{2}{3} \end{array}\right]
\left[\begin{array}{c}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{array}\right]=
\left[\begin{array}{c} 0\\ 1\\ 0\\ 0\\ 0 \end{array}\right].
\end{equation}\]
Therefore,
(4.40)#\[\begin{align}
\left[\begin{array}{c}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{array}\right] & =
\left[\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1\\ -2 & -1 & 0 & 1 & 2\\ 2 & \frac{1}{2} & 0 & \frac{1}{2} & 2\\ -\frac{4}{3} & -\frac{1}{6} & 0 & \frac{1}{6} & \frac{4}{3}\\ \frac{2}{3} & \frac{1}{24} & 0 & \frac{1}{24} & \frac{2}{3} \end{array}\right]^{-1}
\left[\begin{array}{c} 0\\ 1\\ 0\\ 0\\ 0 \end{array}\right]
\notag \\
&= \left[\begin{array}{ccccc} 0 & \frac{1}{12} & -\frac{1}{12} & -\frac{1}{2} & 1\\ 0 & -\frac{2}{3} & \frac{4}{3} & 1 & -4\\ 1 & 0 & -\frac{5}{2} & 0 & 6\\ 0 & \frac{2}{3} & \frac{4}{3} & -1 & -4\\ 0 & -\frac{1}{12} & -\frac{1}{12} & \frac{1}{2} & 1 \end{array}\right]
\left[\begin{array}{c} 0\\ 1\\ 0\\ 0\\ 0 \end{array}\right]
=
\left[\begin{array}{r} \frac{1}{12}\\ -\frac{2}{3}\\ 0\\ \frac{2}{3}\\ -\frac{1}{12} \end{array}\right],
%\stackrel{~or~}{=}\frac{1}{12} \left[\begin{array}{c} 1\\ -8\\ 0\\ 8\\ -1 \end{array}\right]
\end{align}\]
We have,
(4.41)#\[\begin{equation}
f'(x_{i}) = \frac{1}{h}
\left(\frac{1}{12}f_{i-2} - \frac{2}{3}f_{i-1} + \frac{2}{3}f_{i+1} - \frac{1}{12}f_{i+2}\right) + \mathcal{O}(h^4).
\end{equation}\]
Thus,
(4.42)#\[\begin{align}
{f'(x_{i}) \approx \frac{1}{h}\left(\frac{1}{12}f_{i-2} - \frac{2}{3}f_{i-1} + \frac{2}{3}f_{i+1} - \frac{1}{12}f_{i+2}\right).}
\end{align}\]
\subsection{Higher-order approximations of the second derivative of a function} Similarly,
(4.43)#\[\begin{align}\label{2ndder_4thorder_eq01}
h^2 f''(x_i) = \sum_{j = -2}^{2} a_{j+2} f_{i+j}
&= a_{0}f_{i-2} + a_{1}f_{i-1} + a_{2}f_{i} + a_{3}f_{i+1} + a_{4}f_{i+2}
\notag \\ &= a_{0}\left(f_{i} - 2hf_{i}' + 2h^2 f_{i}'' - \frac{4}{3} h^3 f_{i}^{(3)} + \frac{2}{3} h^4 f_{i}^{(4)} \right)
\notag \\ &
+ a_{1}\left( 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)} \right)
\notag \\ &
+ a_{2}f_{i}
\notag \\ &
+ a_{3}\left( 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)} \right)
\notag \\ &
+ a_{4}\left( f_{i} + 2hf_{i}' + 2h^2 f_{i}'' + \frac{4}{3} h^3 f_{i}^{(3)} + \frac{2}{3} h^4 f_{i}^{(4)} \right)+ \mathcal{O}(h^5).
\end{align}\]
Therefore,
(4.44)#\[\begin{align}
h^2 f''(x_i) &= \left(a_{0} + a_{1} + a_{2} + a_{3} +a_{4}\right)f_{i}
\notag \\ &
+ \left(-2a_{0} - a_{1} + 0a_{2} + a_{3} +2a_{4}\right)h\,f_{i}'
\notag \\ &
+ \left(2a_{0} +\frac{1}{2} a_{1} + 0a_{2} + \dfrac{1}{2} a_{3} + 2a_{4}\right)h^2\,f_{i}''
\notag \\ &
+ \left( -\dfrac{4}{3}a_{0} -\dfrac{1}{6} a_{1} + 0a_{2} + \dfrac{1}{6} a_{3} + \dfrac{4}{3}a_{4}\right)h^3\,f_{i}^{(3)}
\notag \\ &
+ \left(\dfrac{2}{3}a_{0} +\dfrac{1}{24} a_{1} + 0a_{2} + \dfrac{1}{24} a_{3} + \dfrac{2}{3}a_{4}\right)h^4\,f_{i}^{(4)}
+ \mathcal{O}(h^5).
\end{align}\]
Since \(h>0\), it follows that
(4.45)#\[\begin{align}
\begin{cases}
a_{0} + a_{1} + a_{2} + a_{3} +a_{4} &=0, \\
-2a_{0} - a_{1} + 0a_{2} + a_{3} +2a_{4} &=0, \\
2a_{0} +\dfrac{1}{2} a_{1} + 0a_{2} + \dfrac{1}{2} a_{3} + 2a_{4} &=1, \\
-\dfrac{4}{3}a_{0} -\dfrac{1}{6} a_{1} + 0a_{2} + \dfrac{1}{6} a_{3} + \dfrac{4}{3}a_{4} &=0, \\
\dfrac{2}{3}a_{0} +\dfrac{1}{24} a_{1} + 0a_{2} + \dfrac{1}{24} a_{3} + \dfrac{2}{3}a_{4} &=0.
\end{cases}
\end{align}\]
The solution to the system can be found as follows,
(4.46)#\[\begin{equation}
\left[\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1\\ -2 & -1 & 0 & 1 & 2\\ 2 & \frac{1}{2} & 0 & \frac{1}{2} & 2\\ -\frac{4}{3} & -\frac{1}{6} & 0 & \frac{1}{6} & \frac{4}{3}\\ \frac{2}{3} & \frac{1}{24} & 0 & \frac{1}{24} & \frac{2}{3} \end{array}\right]
\left[\begin{array}{c}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{array}\right]=
\left[\begin{array}{c} 0\\ 0\\ 1\\ 0\\ 0 \end{array}\right].
\end{equation}\]
Therefore,
(4.47)#\[\begin{align}\label{2ndder_4thorder_eq04}
\left[\begin{array}{c}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{array}\right] & =
\left[\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1\\ -2 & -1 & 0 & 1 & 2\\ 2 & \frac{1}{2} & 0 & \frac{1}{2} & 2\\ -\frac{4}{3} & -\frac{1}{6} & 0 & \frac{1}{6} & \frac{4}{3}\\ \frac{2}{3} & \frac{1}{24} & 0 & \frac{1}{24} & \frac{2}{3} \end{array}\right]^{-1}
\left[\begin{array}{c} 0\\ 0\\ 1\\ 0\\ 0 \end{array}\right]
\notag \\
&= \left[\begin{array}{ccccc} 0 & \frac{1}{12} & -\frac{1}{12} & -\frac{1}{2} & 1\\ 0 & -\frac{2}{3} & \frac{4}{3} & 1 & -4\\ 1 & 0 & -\frac{5}{2} & 0 & 6\\ 0 & \frac{2}{3} & \frac{4}{3} & -1 & -4\\ 0 & -\frac{1}{12} & -\frac{1}{12} & \frac{1}{2} & 1 \end{array}\right]
\left[\begin{array}{c} 0\\ 0\\ 1\\ 0\\ 0 \end{array}\right]
=
\left[\begin{array}{r} -\frac{1}{12}\\ \frac{4}{3}\\ -\frac{5}{2}\\ \frac{4}{3}\\ -\frac{1}{12} \end{array}\right].
\end{align}\]
Using the above coefficients, we have,
(4.48)#\[\begin{equation}\label{2ndder_4thorder_eq02}
f''(x_i) = \frac{1}{h^2}\left(-\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} -\frac{5}{2} f_{i} + \frac{4}{3}f_{i+1} -\frac{1}{12}f_{i+2}\right) + \mathcal{O}(h^4).
\end{equation}\]
Thus,
(4.49)#\[\begin{align}\label{2ndder_4thorder_eq03}
{f''(x_i) \approx \frac{1}{h^2}\left(-\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} -\frac{5}{2} f_{i} + \frac{4}{3}f_{i+1} -\frac{1}{12}f_{i+2}\right).}
\end{align}\]
4.3.2. Higher-order approximations of the second derivative of a function#
Similarly,
(4.50)#\[\begin{align}
h^2 f''(x_i) = \sum_{j = -2}^{2} a_{j+2} f_{i+j}
&= a_{0}f_{i-2} + a_{1}f_{i-1} + a_{2}f_{i} + a_{3}f_{i+1} + a_{4}f_{i+2}
\notag \\ &= a_{0}\left(f_{i} - 2hf_{i}' + 2h^2 f_{i}'' - \frac{4}{3} h^3 f_{i}^{(3)} + \frac{2}{3} h^4 f_{i}^{(4)} \right)
\notag \\ &
+ a_{1}\left( 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)} \right)
\notag \\ &
+ a_{2}f_{i}
\notag \\ &
+ a_{3}\left( 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)} \right)
\notag \\ &
+ a_{4}\left( f_{i} + 2hf_{i}' + 2h^2 f_{i}'' + \frac{4}{3} h^3 f_{i}^{(3)} + \frac{2}{3} h^4 f_{i}^{(4)} \right)+ \mathcal{O}(h^5).
\end{align}\]
Therefore,
(4.51)#\[\begin{align}
h^2 f''(x_i) &= \left(a_{0} + a_{1} + a_{2} + a_{3} +a_{4}\right)f_{i}
\notag \\ &
+ \left(-2a_{0} - a_{1} + 0a_{2} + a_{3} +2a_{4}\right)h\,f_{i}'
\notag \\ &
+ \left(2a_{0} +\frac{1}{2} a_{1} + 0a_{2} + \dfrac{1}{2} a_{3} + 2a_{4}\right)h^2\,f_{i}''
\notag \\ &
+ \left( -\dfrac{4}{3}a_{0} -\dfrac{1}{6} a_{1} + 0a_{2} + \dfrac{1}{6} a_{3} + \dfrac{4}{3}a_{4}\right)h^3\,f_{i}^{(3)}
\notag \\ &
+ \left(\dfrac{2}{3}a_{0} +\dfrac{1}{24} a_{1} + 0a_{2} + \dfrac{1}{24} a_{3} + \dfrac{2}{3}a_{4}\right)h^4\,f_{i}^{(4)}
+ \mathcal{O}(h^5).
\end{align}\]
Since \(h>0\), it follows that
(4.52)#\[\begin{align}
\begin{cases}
a_{0} + a_{1} + a_{2} + a_{3} +a_{4} &=0, \\
-2a_{0} - a_{1} + 0a_{2} + a_{3} +2a_{4} &=0, \\
2a_{0} +\dfrac{1}{2} a_{1} + 0a_{2} + \dfrac{1}{2} a_{3} + 2a_{4} &=1, \\
-\dfrac{4}{3}a_{0} -\dfrac{1}{6} a_{1} + 0a_{2} + \dfrac{1}{6} a_{3} + \dfrac{4}{3}a_{4} &=0, \\
\dfrac{2}{3}a_{0} +\dfrac{1}{24} a_{1} + 0a_{2} + \dfrac{1}{24} a_{3} + \dfrac{2}{3}a_{4} &=0.
\end{cases}
\end{align}\]
The solution to the system can be found as follows,
(4.53)#\[\begin{equation}
\left[\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1\\ -2 & -1 & 0 & 1 & 2\\ 2 & \frac{1}{2} & 0 & \frac{1}{2} & 2\\ -\frac{4}{3} & -\frac{1}{6} & 0 & \frac{1}{6} & \frac{4}{3}\\ \frac{2}{3} & \frac{1}{24} & 0 & \frac{1}{24} & \frac{2}{3} \end{array}\right]
\left[\begin{array}{c}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{array}\right]=
\left[\begin{array}{c} 0\\ 0\\ 1\\ 0\\ 0 \end{array}\right].
\end{equation}\]
Therefore,
(4.54)#\[\begin{align}
\left[\begin{array}{c}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\end{array}\right] & =
\left[\begin{array}{ccccc} 1 & 1 & 1 & 1 & 1\\ -2 & -1 & 0 & 1 & 2\\ 2 & \frac{1}{2} & 0 & \frac{1}{2} & 2\\ -\frac{4}{3} & -\frac{1}{6} & 0 & \frac{1}{6} & \frac{4}{3}\\ \frac{2}{3} & \frac{1}{24} & 0 & \frac{1}{24} & \frac{2}{3} \end{array}\right]^{-1}
\left[\begin{array}{c} 0\\ 0\\ 1\\ 0\\ 0 \end{array}\right]
\notag \\
&= \left[\begin{array}{ccccc} 0 & \frac{1}{12} & -\frac{1}{12} & -\frac{1}{2} & 1\\ 0 & -\frac{2}{3} & \frac{4}{3} & 1 & -4\\ 1 & 0 & -\frac{5}{2} & 0 & 6\\ 0 & \frac{2}{3} & \frac{4}{3} & -1 & -4\\ 0 & -\frac{1}{12} & -\frac{1}{12} & \frac{1}{2} & 1 \end{array}\right]
\left[\begin{array}{c} 0\\ 0\\ 1\\ 0\\ 0 \end{array}\right] =
\left[\begin{array}{r} -\frac{1}{12}\\ \frac{4}{3}\\ -\frac{5}{2}\\ \frac{4}{3}\\ -\frac{1}{12} \end{array}\right].
\end{align}\]
Using the above coefficients, we have,
(4.55)#\[\begin{equation}
f''(x_i) = \frac{1}{h^2}\left(-\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} -\frac{5}{2} f_{i} + \frac{4}{3}f_{i+1} -\frac{1}{12}f_{i+2}\right) + \mathcal{O}(h^4).
\end{equation}\]
Thus,
(4.56)#\[\begin{align}
{f''(x_i) \approx \frac{1}{h^2}\left(-\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} -\frac{5}{2} f_{i} + \frac{4}{3}f_{i+1} -\frac{1}{12}f_{i+2}\right).}
\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\). Then,
a. Approximate \(f'(x)\) and \(f''(x)\) using the above fourth-order methods at
\[\begin{equation*}
\{0.02,~0.03,~0.04, \ldots,~0.97,~0.98\}.
\end{equation*}\]
b. Demonstrate their order of convergence numerically.
Solution:
a.
\[\begin{align*}
f'(x_i) & \approx \frac{1}{h}\left(\frac{1}{12}f_{i-2} - \frac{2}{3}f_{i-1} - \frac{2}{3}f_{i+1} + \frac{1}{12}f_{i+2}\right),\\
f''(x_i) & \approx \frac{1}{h^2}\left(-\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} -\frac{5}{2} f_{i} + \frac{4}{3}f_{i+1} -\frac{1}{12}f_{i+2}\right).
\end{align*}\]
For \hl{\(x_{2} = 0.02\)}, since \(f(x_i+ jh)\) and \(h=0.01\), we have,
\[\begin{align*}
f'(0.02) & \approx \frac{1}{h}\left( \frac{1}{12}f(0.02 -2h) - \frac{2}{3}f(0.02 -h) - \frac{2}{3}f(0.02 + h) + \frac{1}{12}f(0.02 + 2h) \right)
\\ & = \frac{1}{0.01}\left( \frac{1}{12}f(0) - \frac{2}{3}f(0.01) - \frac{2}{3}f(0.03) + \frac{1}{12}f(0.04) \right) = 3.141815.
\end{align*}\]
Similarly,
\[\begin{align*}
f'(0.03) & \approx \frac{1}{0.01}\left( \frac{1}{12}f(0.01) - \frac{2}{3}f(0.02) - \frac{2}{3}f(0.04) + \frac{1}{12}f(0.05) \right) = 3.214100,\\
f'(0.04) & \approx \frac{1}{0.01}\left( \frac{1}{12}f(0.02) - \frac{2}{3}f(0.03) - \frac{2}{3}f(0.05) + \frac{1}{12}f(0.06) \right) = 3.287319,\\
\vdots &\\
f'(0.97) & \approx \frac{1}{0.01}\left( \frac{1}{12}f(0.95) - \frac{2}{3}f(0.96) - \frac{2}{3}f(0.98) + \frac{1}{12}f(0.99) \right) = 16.415137,\\
f'(0.98) & \approx \frac{1}{0.01}\left( \frac{1}{12}f(0.96) - \frac{2}{3}f(0.97) - \frac{2}{3}f(0.99) + \frac{1}{12}f(1.00) \right) = 16.657367.
\end{align*}\]
import numpy as np
import pandas as pd
f = lambda x: 3*x*np.exp(x) - np.cos(x)
f1=lambda x: 3*np.exp(x) + 3*x*np.exp(x) + np.sin(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)
\[\begin{align*}
E_{h}(\text{Fourth-order approximation of the first derivative}) &= \max_{i}\left|
f'(x) - \frac{1}{h}\left(\frac{1}{12}f_{i-2} - \frac{2}{3}f_{i-1} - \frac{2}{3}f_{i+1} + \frac{1}{12}f_{i+2}\right)
\right|,\\
E_{h}(\text{Fourth-order approximation of the second derivative})&= \max_{i}\left|
f''(x) - \frac{1}{h^2}\left(-\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} -\frac{5}{2} f_{i} + \frac{4}{3}f_{i+1} -\frac{1}{12}f_{i+2}\right)
\right|.
\end{align*}\]
from IPython.display import display, Latex
Ind = np.arange(2, len(xn)-2)
# Higher-order approximations of the first derivative of a function
F1App = ((1/12)*yn[Ind-2] - (2/3)* yn[Ind-1] + (2/3)* yn[Ind+1] - (1/12)*yn[Ind+2]) / h
# f'(x) Exact
F1_Exact = f1(xn[Ind])
# Higher-order approximations of the second derivative of a function
F2App = (-(1/12)*yn[Ind-2] + (4/3)* yn[Ind-1] -(5/2)* yn[Ind] + (4/3)* yn[Ind+1] - (1/12)*yn[Ind+2]) / (h**2)
# f''(x) Exact
F2_Exact = f2(xn[Ind])
display(Latex('''\\max_{i}\\left|
f'(x) - \\frac{1}{h}\\left(\\frac{1}{12}f_{i-2} - \\frac{2}{3}f_{i-1} - \\frac{2}{3}f_{i+1} + \\frac{1}{12}f_{i+2}\\right)
\\right| = %.4e''' % max(abs(F1App - F1_Exact))))
display(Latex('''\\max_{i}\\left|
f''(x) - \\frac{1}{h^2}\\left(-\\frac{1}{12}f_{i-2} + \\frac{4}{3}f_{i-1} -\\frac{5}{2} f_{i} + \\frac{4}{3}f_{i+1} -\\frac{1}{12}f_{i+2}\\right)
\\right| = %.4e''' % max(abs(F2App - F2_Exact))))
\[\max_{i}\left|
f'(x) - \frac{1}{h}\left(\frac{1}{12}f_{i-2} - \frac{2}{3}f_{i-1} - \frac{2}{3}f_{i+1} + \frac{1}{12}f_{i+2}\right)
\right| = 1.6211e-08\]
\[\max_{i}\left|
f''(x) - \frac{1}{h^2}\left(-\frac{1}{12}f_{i-2} + \frac{4}{3}f_{i-1} -\frac{5}{2} f_{i} + \frac{4}{3}f_{i+1} -\frac{1}{12}f_{i+2}\right)
\right| = 6.2761e-09\]
The order of convergence numerically:
h = [2**(-i) for i in range(3, 8)]
Cols = ['h', 'Max Error (First Derivative)', 'Max Error (Second Derivative)']
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)
# Innter points
Ind = np.arange(2, len(xn)-2)
# Exact
F1Exact = f1(xn[Ind])
F2Exact = f2(xn[Ind])
# A Fourth-order Approximation of the first derivative
F1App = ((1/12)*yn[Ind-2] - (2/3)* yn[Ind-1] + (2/3)* yn[Ind+1] - (1/12)*yn[Ind+2]) / h[i]
Table[Cols[1]][i] = np.max(np.abs(F1Exact - F1App))
# A Fourth-order Approximation of the Second derivative
F2App = (-(1/12)*yn[Ind-2] + (4/3)* yn[Ind-1] -(5/2)* yn[Ind] + (4/3)* yn[Ind+1] - (1/12)*yn[Ind+2]) / (h[i]**2)
Table[Cols[2]][i] = np.max(np.abs(F2Exact - F2App))
del Ind, F1Exact, F1App, F2Exact, F2App
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 | Max Error (First Derivative) | Max Error (Second Derivative) | |
|---|---|---|---|---|
| 0 | 0.125000 | 8.0000e+00 | 2.6196e-04 | 1.0311e-04 |
| 1 | 0.062500 | 1.6000e+01 | 2.0369e-05 | 7.9286e-06 |
| 2 | 0.031250 | 3.2000e+01 | 1.4193e-06 | 5.4997e-07 |
| 3 | 0.015625 | 6.4000e+01 | 9.3660e-08 | 3.6226e-08 |
| 4 | 0.007812 | 1.2800e+02 | 6.0149e-09 | 2.3790e-09 |
# 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['Max Error (First Derivative)'].values,
Table['Max Error (Second Derivative)'].values],
labels = ['First Derivative', 'Second Derivative'],
title = 'Order of accuracy: First and Second Derivatives',
legend_orientation = 'horizontal', ylim = [3.5, 4.1])