Processing math: 100%

Chapter 1
HWs

1.1 HW1

PDF (letter size)

PDF (legal size)

1.1.1 Problem 1

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 )$.

1.1.2 Problem 2

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). Hence¹

$$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

f[x_] := (8 Sin[x] - Sin[2 x])/6
Plot[{x, f[x]}, {x, 0, Pi}, Frame -> True,
FrameTicks -> {{Automatic, None}, {Range[0, Pi, Pi/4], None}},
FrameLabel -> {{Text@Style["effective k dx"],
None}, {Text@Style["Actual k dx"],
Text@Style[
Column[{"Effective wave numbers",
"First derivative Imaginary component"}, Alignment -> Center],
Bold]}}, BaseStyle -> 14,
PlotLegends -> {"Exact", "3rd order"}, GridLines -> Automatic,
GridLinesStyle -> Directive[LightGray, Dashed],
PlotStyle -> {Red, Blue}]

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

f[x_] := (Cos[2 x] - 4 Cos[x] + 3)/6
Plot[{0, f[x]}, {x, 0, Pi}, Frame -> True,
FrameTicks -> {{Automatic, None}, {Range[0, Pi, Pi/4], None}},
FrameLabel -> {{Text@Style["effective k dx"],
None}, {Text@Style["Actual k dx"],
Text@Style[
Column[{"Effective wave numbers",
"First derivative Real component"}, Alignment -> Center],
Bold]}}, BaseStyle -> 14, PlotLegends -> {"Exact", "3rd order"},
GridLines -> Automatic,
GridLinesStyle -> Directive[LightGray, Dashed],
PlotStyle -> {Red, Blue}]

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.

real[x_] := (Cos[2 x] - 4 Cos[x] + 3)/6
im[x_] := (8 Sin[x] - Sin[2 x])/6
Plot[{real[x], im[x]}, {x, 0, Pi}, Frame -> True,
FrameTicks -> {{Automatic, None}, {Range[0, Pi, Pi/4], None}},
FrameLabel -> {{Text@Style["effective k dx"],
None}, {Text@Style["Actual k dx"],
Text@Style[
Column[{"Effective wave numbers",
"First derivative Real vs. Imaginary components"},
Alignment -> Center], Bold]}}, BaseStyle -> 14,
PlotLegends -> {"Real 3rd order", "Imaginary 3rd order"},
GridLines -> Automatic,
GridLinesStyle -> Directive[LightGray, Dashed],
PlotStyle -> {Red, Blue}]

1.1.3 Problem 3

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.