Let \begin{equation} \frac{du}{dx}\approx \left . \frac{\delta u}{\delta x}\right \vert _{i}=au_{i-2}+bu_{i-1}+cu_{i}+du_{i+1} \tag{1} \end{equation} We now set up the Taylor table as explained in the lecture notes using \(h\) in place of \(dx\) for the spatial grid spacing in order to simplify the notation. Since we want to find \(4\) unknowns (\(a,b,c,d\)), then we need at least 4 columns. But we generate 5 in order to check for the order of the error using the last column. Therefore, the Taylor table with 5 columns is
| \(u_{i}\) | \(\left . \frac{\partial u}{\partial x}\right \vert _{i}\) | \(\left . \frac{\partial ^{2}u}{\partial x^{2}}\right \vert _{i}\) | \(\left . \frac{\partial ^{3}u}{\partial x^{3}}\right \vert _{i}\) | \(\left . \frac{\partial ^{4}u}{\partial x^{4}}\right \vert _{i}\) | |
| \(u_{i-2}\) | \(1\) | \(-2h\) | \(\left ( -2h\right ) ^{2}\frac{1}{2!}\) | \(\left ( -2h\right ) ^{3}\frac{1}{3!}\) | \(\left ( -2h\right ) ^{4}\frac{1}{4!}\) |
| \(u_{i-1}\) | \(1\) | \(-h\) | \(\left ( -h\right ) ^{2}\frac{1}{2!}\) | \(\left ( -h\right ) ^{3}\frac{1}{3!}\) | \(\left ( -h\right ) ^{4}\frac{1}{4!}\) |
| \(u_{i}\) | \(1\) | \(0\) | \(0\) | \(0\) | \(0\) |
| \(u_{i+1}\) | \(1\) | \(h\) | \(h^{2}\frac{1}{2!}\) | \(h^{3}\frac{1}{3!}\) | \(h^{4}\frac{1}{4!}\) |
We now add the coefficients \(a,b,c,\) and \(d\) to obtain
| \(u_{i}\) | \(\left . \frac{\partial u}{\partial x}\right \vert _{i}\) | \(\left . \frac{\partial ^{2}u}{\partial x^{2}}\right \vert _{i}\) | \(\left . \frac{\partial ^{3}u}{\partial x^{3}}\right \vert _{i}\) | \(\left . \frac{\partial ^{4}u}{\partial x^{4}}\right \vert _{i}\) | |
| \(au_{i-2}\) | \(a\) | \(a\left ( -2h\right ) \) | \(a\left ( -2h\right ) ^{2}\frac{1}{2!}\) | \(a\left ( -2h\right ) ^{3}\frac{1}{3!}\) | \(a\left ( -2h\right ) ^{4}\frac{1}{4!}\) |
| \(bu_{i-1}\) | \(b\) | \(b\left ( -h\right ) \) | \(b\left ( -h\right ) ^{2}\frac{1}{2!}\) | \(b\left ( -h\right ) ^{3}\frac{1}{3!}\) | \(b\left ( -h\right ) ^{4}\frac{1}{4!}\) |
| \(cu_{i}\) | \(c\) | \(0\) | \(0\) | \(0\) | \(0\) |
| \(du_{i+1}\) | \(d\) | \(d\left ( h\right ) \) | \(d\left ( h\right ) ^{2}\frac{1}{2!}\) | \(d\left ( h\right ) ^{3}\frac{1}{3!}\) | \(d\left ( h\right ) ^{4}\frac{1}{4!}\) |
Expanding and summing each column gives
| \(u_{i}\) | \(\left . \frac{\partial u}{\partial x}\right \vert _{i}\) | \(\left . \frac{\partial ^{2}u}{\partial x^{2}}\right \vert _{i}\) | \(\left . \frac{\partial ^{3}u}{\partial x^{3}}\right \vert _{i}\) | \(\left . \frac{\partial ^{4}u}{\partial x^{4}}\right \vert _{i}\) | |
| \(au_{i-2}\) | \(a\) | \(a\left ( -2h\right ) \) | \(a\left ( -2h\right ) ^{2}\frac{1}{2!}\) | \(a\left ( -2h\right ) ^{3}\frac{1}{3!}\) | \(a\left ( -2h\right ) ^{4}\frac{1}{4!}\) |
| \(bu_{i-1}\) | \(b\) | \(b\left ( -h\right ) \) | \(b\left ( -h\right ) ^{2}\frac{1}{2!}\) | \(b\left ( -h\right ) ^{3}\frac{1}{3!}\) | \(b\left ( -h\right ) ^{4}\frac{1}{4!}\) |
| \(cu_{i}\) | \(c\) | \(0\) | \(0\) | \(0\) | \(0\) |
| \(du_{i+1}\) | \(d\) | \(d\left ( h\right ) \) | \(d\left ( h\right ) ^{2}\frac{1}{2!}\) | \(d\left ( h\right ) ^{3}\frac{1}{3!}\) | \(d\left ( h\right ) ^{4}\frac{1}{4!}\) |
| \(\sum \) | \(a+b+c+d\) | \(\left ( -2a-b+d\right ) h\) | \(\left ( 2a+\frac{b}{2}+\frac{d}{2}\right ) h^{2}\) | \(\left ( -\frac{8}{6}a-\frac{b}{6}+\frac{d}{6}\right ) h^{3}\) | \(\left ( \frac{16}{24}a+\frac{b}{24}+\frac{d}{24}\right ) h^{4}\) |
| \(0\) | \(1\) | \(0\) | \(0\) | check if zero | |
Since first derivative approximation is sought, we want the \(\frac{\partial u}{\partial x}\) column to sum to one, and the other columns to sum to zero. This gives four equations to solve for \(a,b,c\) and \(d\)\begin{align*} a+b+c+d & =0\\ \left ( -2a-b+d\right ) h & =1\\ \left ( 2a+\frac{b}{2}+\frac{d}{2}\right ) h^{2} & =0\\ \left ( -\frac{8}{6}a-\frac{b}{6}+\frac{d}{6}\right ) h^{3} & =0 \end{align*}
Since \(h\neq 0\) these reduce to\begin{align*} a+b+c+d & =0\\ -2a-b+d & =\frac{1}{h}\\ 2a+\frac{b}{2}+\frac{d}{2} & =0\\ -\frac{8}{6}a-\frac{b}{6}+\frac{d}{6} & =0 \end{align*}
Solving gives \(a=\frac{1}{6h},b=-\frac{1}{h},c=\frac{1}{2h},d=\frac{1}{3h}\). Therefore (1) becomes\begin{align*} \left . \frac{du}{dx}\right \vert _{x_{i}} & \approx \left . \frac{\delta u}{\delta x}\right \vert _{i}=au_{i-2}+bu_{i-1}+cu_{i}+du_{i+1}\\ & =\frac{\frac{1}{6}u_{i-2}-u_{i-1}+\frac{1}{2}u_{i}+\frac{1}{3}u_{i+1}}{h}\\ & =\frac{u_{i-2}-6u_{i-1}+3u_{i}+2u_{i+1}}{6h} \end{align*}
To determine the truncation error the last column in the Taylor table above is checked if it sums to non-zero. If the sum turns out to be zero, the next column after that must then be checked.\begin{align*} \left ( \frac{16}{24}a+\frac{b}{24}+\frac{d}{24}\right ) h^{4} & =\left ( \frac{16}{24}\frac{1}{6h}-\frac{1}{24h}+\frac{1}{3\left ( 24\right ) h}\right ) h^{4}\\ & =\left ( \frac{16}{24}\frac{1}{6}-\frac{1}{24}+\frac{1}{3\left ( 24\right ) }\right ) h^{3}\\ & =\frac{1}{12}h^{3} \end{align*}
Since the sum is not zero, there is no need to check any more columns and the truncation error is verified to be third order \(O\left ( h^{3}\right ) \).
Using result from problem 1\begin{equation} \left . \frac{\delta u}{\delta x}\right \vert _{i}=\frac{u_{i-2}-6u_{i-1}+3u_{i}+2u_{i+1}}{6h} \tag{1} \end{equation} Using \[ u\left ( x\right ) =\sum _{k}\hat{u}_{k}e^{jkx}\] Where \(\hat{u}_{k}\) are the Fourier coefficients, which are functions of \(k\), and are complex numbers in general . Looking at one mode only (one specific \(k\)), then we let \(k\) run over its range, where \(k\) is called the wave number which is related to the wave length \(\lambda \) by\[ k=\frac{2\pi }{\lambda }\] \(j\) above is \(\sqrt{-1}\)(We could also have used \(\hat{\imath }\) for \(\sqrt{-1}\) but it looked very close to the index \(i\) and can be confusing). Hence1
\[ u\left ( x\right ) =\hat{u}_{k}e^{jkx}\] Equation (1) now can be written as\begin{align} \frac{\partial u}{\partial x} & =\frac{\partial }{\partial x}\hat{u}_{k}e^{jkx}\nonumber \\ & =\left ( jk\right ) \hat{u}_{k}e^{jkx}\nonumber \\ & =\left ( jk\right ) u\left ( x\right ) \tag{2} \end{align}
For finite difference the above can be written as
\[ \left . \frac{\delta u}{\delta x}\right \vert _{i}=\left ( jk\right ) _{eff}u_{i}\]
And the goal is to determine \(\left ( jk\right ) _{eff}\) using (1) above and compare it to the actual \(\left ( jk\right ) \) from (2). From (1) we obtain for the RHS\begin{align*} \left ( jk\right ) _{eff}u_{i} & =\frac{\hat{u}_{k}e^{jk\left ( x_{i}-2h\right ) }-6\hat{u}_{k}e^{jk\left ( x_{i}-h\right ) }+3\hat{u}_{k}e^{jkx_{i}}+2\hat{u}_{k}e^{jk\left ( x_{i}+h\right ) }}{6h}\\ \left ( jk\right ) _{eff}u_{i} & =\left ( \frac{e^{-2jkh}-6e^{-jkh}+3+2e^{jkh}}{6h}\right ) \hat{u}_{k}e^{jkx_{i}}\\ \left ( jk\right ) _{eff}u_{i} & =\overset{\text{effective wave number}}{\overbrace{\frac{e^{-2jkh}-6e^{-jkh}+3+2e^{jkh}}{6h}}}u_{i} \end{align*}
Therefore the effective wave number\(\left ( jk\right ) _{eff}\) is\begin{align*} \left ( jk\right ) _{eff} & =\frac{e^{-2jkh}-6e^{-jkh}+3+2e^{jkh}}{6h}\\ & =\frac{\left ( \cos 2kh-j\sin 2kh\right ) -6\left ( \cos kh-j\sin kh\right ) +3+2\left ( \cos kh+j\sin kh\right ) }{6h}\\ & =\frac{j}{6h}\left ( -\sin 2kh+6\sin kh+2\sin kh\right ) +\frac{1}{6h}\left ( \cos 2kh-6\cos kh+3+2\cos kh\right ) \end{align*}
Therefore\[ \left ( jk\right ) _{eff}=j\overset{\text{complex part}}{\overbrace{\left ( \frac{8\sin kh-\sin 2kh}{6h}\right ) }}+\overset{\text{real part}}{\overbrace{\frac{1}{6h}\left ( \cos 2kh-4\cos kh+3\right ) }}\] We see that \(\left ( jk\right ) _{eff}\) has both a complex part and a real part. But the exact wave number \(\left ( jk\right ) \) is only complex. This is the first major difference we see. Now we will plot the real and the imaginary parts of \(\left ( jk\right ) _{eff}.\) The complex part is\[ \left ( jk\right ) _{eff_{\operatorname{complex}}}=\frac{8\sin kh-\sin 2kh}{6}\] And the second is the real part\[ \left ( jk\right ) _{eff_{\operatorname{real}}}=\frac{\cos 2kh-4\cos kh+3}{6}\] We now use \(x\) for \(kh\) as the argument to simplify the notation and plot it\[ k_{eff_{\operatorname{complex}}}\left ( x\right ) =\frac{8\sin x-\sin 2x}{6}\] And the real part is\[ k_{eff_{\operatorname{real}}}\left ( x\right ) =\frac{\cos 2x-4\cos x+3}{6}\] The plots of the imaginary part is given below
Discussion: We see from the above that the imaginary part of the effective wave number is accurate and close to the exact value for small wave numbers. After about \(kh\approx \frac{\pi }{3}\), then it is no longer accurate. Smaller \(k\) implies larger wave length \(\lambda \) which in turn puts a limits of the grid size \(h.\)
The real part plot is below
Discussion: The exact value is zero for all wave numbers, since we know from the above, that the exact effective \(k\) has only complex part and no real part. but the effective \(k\) is only as accurate and close to zero for much smaller wave numbers. After about \(kh\approx \frac{\pi }{4}\) it is no longer accurate. Having a real part in the effective wave number, implies the finite difference scheme will introduce damping effect in the result.
We need to derive approximation for \(\left . \frac{du}{dx}\right \vert _{x_{i}}\approx \left . \frac{\delta u}{\delta x}\right \vert _{i}=\frac{u_{i+1/2}\left ( x\right ) -u_{i-1/2}\left ( x\right ) }{h}\) using 3 points Lagrangian interpolation. There are 4 points needed. The following diagram shows the cell structure used
When interpolating \(u_{i+1/2}\left ( x\right ) \), the following 3 points are used
When interpolating for \(u_{i-1/2}\left ( x\right ) \), the following 3 points are used
Therefore
\begin{align*} u_{i+1/2}\left ( x\right ) & =u_{i-1}\left ( \cdot \right ) +u_{i}\left ( \cdot \right ) +u_{i+1}\left ( \cdot \right ) \\ & =u_{i-1}\frac{\left ( x-x_{i}\right ) \left ( x-x_{i+1}\right ) }{\left ( x_{i-1}-x_{i}\right ) \left ( x_{i-1}-x_{i+1}\right ) }+u_{i}\frac{\left ( x-x_{i-1}\right ) \left ( x-x_{i+1}\right ) }{\left ( x_{i}-x_{i-1}\right ) \left ( x_{i}-x_{i+1}\right ) }+u_{i+1}\frac{\left ( x-x_{i-1}\right ) \left ( x-x_{i}\right ) }{\left ( x_{i+1}-x_{i-1}\right ) \left ( x_{i+1}-x_{i}\right ) } \end{align*}
When \(x\) is midpoint between \(x_{i+1}\) and \(x_{i}\), then the above reduces to (where \(h=dx\)) which is the grid size between each point:\begin{align*} u_{i+1/2}\left ( x\right ) & =u_{i-1}\frac{\left ( \frac{h}{2}\right ) \left ( \frac{-h}{2}\right ) }{\left ( -h\right ) \left ( -2h\right ) }+u_{i}\frac{\left ( \frac{3}{2}h\right ) \left ( \frac{-h}{2}\right ) }{\left ( h\right ) \left ( -h\right ) }+u_{i+1}\frac{\left ( \frac{3}{2}h\right ) \left ( \frac{h}{2}\right ) }{\left ( 2h\right ) \left ( h\right ) }\\ & =-\frac{1}{8}u_{i-1}+\frac{3}{4}u_{i}+\frac{3}{8}u_{i+1} \end{align*}
And\begin{align*} u_{i-1/2}\left ( x\right ) & =u_{i-2}\left ( \cdot \right ) +u_{i-1}\left ( \cdot \right ) +u_{i}\left ( \cdot \right ) \\ & =u_{i-2}\frac{\left ( x-x_{i-1}\right ) \left ( x-x_{i}\right ) }{\left ( x_{i-2}-x_{i-1}\right ) \left ( x_{i-2}-x_{i}\right ) }+u_{i-1}\frac{\left ( x-x_{i-2}\right ) \left ( x-x_{i}\right ) }{\left ( x_{i-1}-x_{i-2}\right ) \left ( x_{i-1}-x_{i}\right ) }+u_{i}\frac{\left ( x-x_{i-2}\right ) \left ( x-x_{i-1}\right ) }{\left ( x_{i+1}-x_{i-2}\right ) \left ( x_{i}-x_{i-1}\right ) } \end{align*}
When \(x\) is midpoint between \(x_{i}\) and \(x_{i-1}\), then the above reduces to (where \(h=dx\)) which is the grid size between each point:\begin{align*} u_{i-1/2}\left ( x\right ) & =u_{i-2}\frac{\left ( \frac{h}{2}\right ) \left ( -\frac{h}{2}\right ) }{\left ( -h\right ) \left ( -2h\right ) }+u_{i-1}\frac{\left ( \frac{3}{2}h\right ) \left ( \frac{-h}{2}\right ) }{\left ( h\right ) \left ( -h\right ) }+u_{i}\frac{\left ( \frac{3}{2}h\right ) \left ( \frac{h}{2}\right ) }{\left ( 2h\right ) \left ( h\right ) }\\ & =-\frac{1}{8}u_{i-2}+\frac{3}{4}u_{i-1}+\frac{3}{8}u_{i} \end{align*}
Therefore\begin{align*} \left . \frac{\delta u}{\delta x}\right \vert _{i} & =\frac{u_{i+1/2}\left ( x\right ) -u_{i-1/2}\left ( x\right ) }{dx}\\ & =\frac{\left ( -\frac{1}{8}u_{i-1}+\frac{3}{4}u_{i}+\frac{3}{8}u_{i+1}\right ) -\left ( -\frac{1}{8}u_{i-2}+\frac{3}{4}u_{i-1}+\frac{3}{8}u_{i}\right ) }{h}\\ & =\frac{1}{8}\frac{3u_{i}-7u_{i-1}+3u_{i+1}+u_{i-2}}{h} \end{align*}
To determine the Taylor series accuracy, we expand the RHS around \(x_{i}\)\begin{align*} \Delta & =\frac{1}{8h}\left ( 3u_{i}-7u_{i-1}+3u_{i+1}+u_{i-2}\right ) \\ & \approx \frac{1}{8h}\left [ 3u_{i}-7\left ( u_{i}-h\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( \left ( -h\right ) ^{2}\right ) \right ) +3\left ( u_{i}+h\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( h^{2}\right ) \right ) +\left ( u_{i}-2h\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( \left ( -2h\right ) ^{2}\right ) \right ) \right ] \\ & =\frac{1}{8h}\left [ 3u_{i}-7u_{i}+7h\left . \frac{\delta u}{\delta x}\right \vert _{i}+7O\left ( h^{2}\right ) +3u_{i}+3h\left . \frac{\delta u}{\delta x}\right \vert _{i}+3O\left ( h^{2}\right ) +u_{i}-2h\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( 4h^{2}\right ) \right ] \\ & =\frac{1}{8h}\left [ 7h\left . \frac{\delta u}{\delta x}\right \vert _{i}+7O\left ( h^{2}\right ) +3h\left . \frac{\delta u}{\delta x}\right \vert _{i}+3O\left ( h^{2}\right ) -2h\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( 4h^{2}\right ) \right ] \\ & =\frac{1}{8h}\left ( 7h\left . \frac{\delta u}{\delta x}\right \vert _{i}+3h\left . \frac{\delta u}{\delta x}\right \vert _{i}-2h\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( h^{2}\right ) \right ) \\ & =\frac{1}{8}\left ( 8\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( h\right ) \right ) \\ & =\left . \frac{\delta u}{\delta x}\right \vert _{i}+O\left ( h\right ) \end{align*}
Therefore this is first order accurate.