This post will discuss solution methods for hyperbolic PDEs. Stability analyses via the method of modified equations are provided alongside code examples.
This may look awful if you are viewing it in dark mode since the notebook container currently forces light mode. This is due to issues with the equations appearing as black text on a dark background.
import matplotlib.pyplot as plt
import numpy as np
Problem Statement¶
Consider the following 1-D hyperbolic PDE: \begin{split} \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 \end{split}
The following initial condition is used, along with periodic boundary conditions.
C = 1
L = 1
NX = 161
dx = L / (NX - 1)
xarr = np.linspace(0, 1, NX)
T = 1.0
NT = 201
Tarr = np.linspace(0, 1, NT)
dt = Tarr[1]
CFL = C * dt / dx
print(f"Simulation CFL number is {CFL:.2f}")
u0 = np.zeros_like(xarr)
u0[NX//16:NX//8] = [16 * dx * i for i in range(NX//16)]
u0[NX//8:NX//4] = C
u0[NX//4:NX//4+NX//16] = [1 - 16 * dx * i for i in range(NX//16)]
plt.figure(figsize=(7, 2))
plt.plot(xarr, u0)
plt.xlabel("x")
plt.ylabel("u(x)")
plt.title("Initial Condition")
Simulation CFL number is 0.80
Text(0.5, 1.0, 'Initial Condition')
1: FTCS¶
For a finite difference method that is Forward in Time and Central in Space (FTCS), the PDE is discretized as: \begin{split} \frac{u_l^{n+1} - u_l^n}{\Delta t} + c \frac{u_{l+1}^n - u_{l-1}^n}{2\Delta x} = 0 \end{split} which results in the following: \begin{split} u_l^{n+1} = u_l^n + \Delta t \bigg( -c \frac{u_{l+1}^n - u_{l-1}^n}{2\Delta x} \bigg) \end{split}
Stability Analysis¶
Assuming a solution of form $u_l^n = A_n e^{ikx_l}$: \begin{split} A_{n+1}e^{ikx_l} = A_ne^{ikx_l} + \Delta t \bigg( -c \frac{A_n e^{ik(x_l + \Delta x)} - A_n e^{ik(x_l - \Delta x)}}{2\Delta x} \bigg) \end{split} Divide out $e^{ikx_l}$: \begin{split} A_{n+1} = A_n - \frac{c \Delta t}{2\Delta x} A_n \big( e^{ik\Delta x} - e^{-ik\Delta x} \big) \end{split}
Using the follow expansions (from Euler's identity): \begin{split} e^{ik\Delta x} &= \cos(k\Delta x) + i \sin(k\Delta x) \\ e^{-ik\Delta x} &= \cos(k\Delta x) - i \sin(k\Delta x) \end{split}
The following is obtained: \begin{split} \frac{A_{n+1}}{A_n} = 1 - \frac{c\Delta t}{\Delta x}i\sin(k\Delta x) \end{split}
The magnitude of the amplification factor is $1 + (\frac{c \Delta t}{\Delta x})^2\sin^2(k\Delta x)$, which is always greater than 1. So, this method is unconditionally unstable.
Modified Equation Analysis¶
Since this method is first order in time and second order in space, the following expansions of each term will be used: \begin{split} u_l^n &= u \\ u_l^{n+1} &= u + \Delta t u_t + \frac{\Delta t^2}{2}u_{tt} + \cdots \\ u_{l+1}^n &= u + \Delta x u_x + \frac{\Delta x^2}{2}u_{xx} + \frac{\Delta x^3}{6}u_{xxx} + \cdots \\ u_{l-1}^n &= u - \Delta x u_x + \frac{\Delta x^2}{2}u_{xx} - \frac{\Delta x^3}{6}u_{xxx} + \cdots \\ \end{split}
Plugging these into the discretized PDE: \begin{split} \frac{\Delta t u_t + \frac{\Delta t^2}{2}u_{tt}}{\Delta t} + \frac{c}{2\Delta x}\bigg(2\Delta x u_x + \frac{\Delta x3}{3}u_{xxx} \bigg) &= 0 \end{split}
\begin{split}\downarrow\end{split}
\begin{split} u_t + \frac{\Delta t}{2}u_{tt} + c\bigg(u_x + \frac{\Delta x^2}{6}u_{xxx}\bigg)=0 \end{split}
To get rid of the double time derivative, use the following: $u_t=-cu_x$, so $u_{tt} = -cu_{xt} = -c(u_t)_x = -c(-cu_x)_x = c^2u_{xx}$. So: \begin{split} u_t + cu_x = -\frac{c^2\Delta t}{2}u_{xx} - \frac{c \Delta x^2}{6}u_{xxx} \end{split}
The diffusion coefficient of $-\frac{c^2\Delta t}{2}$ indicates the equation is actually anti-diffusive, which explains why this method is unconditionally unstable.
This is demonstrated by the simulation below - the solution starts to show nonphysical behaviors after two timesteps (0.02 seconds) and by six timesteps (0.06 second) the simulation is complete nonsense, containing wild oscillations. There is no point in displaying anything past this timestep - these oscillations grow unbounded.
solution_FTCS = np.zeros((NX, NT))
solution_FTCS[:, 0] = u0
for tdx in range(NT - 1):
for xdx in range(NX):
lp1 = (xdx + 1) % NX
lm1 = (xdx - 1) % NX
solution_FTCS[xdx, tdx+1] = (
solution_FTCS[xdx, tdx]
+ dt * (-C * (solution_FTCS[lp1, tdx] - solution_FTCS[lm1, tdx])/(2*dx))
)
plt.figure(figsize=(7,2))
for tdx in [0, 5, 10, 15]:
plt.plot(xarr, solution_FTCS[:, tdx], label=f"{Tarr[tdx]:.2f} s")
plt.legend()
plt.title("FTCS")
Text(0.5, 1.0, 'FTCS')
2: Implicit in Time, Central in Space¶
For a finite difference method that is Backwards in Time and Central in Space (BTCS?), the PDE is discretized as: \begin{split} \frac{u_l^n - u_l^{n-1}}{\Delta t} + c \frac{u_{l+1}^n - u_{l-1}^n}{2\Delta x} = 0 \end{split}
Stability Analysis¶
Assuming a solution of form $u_l^n = A_n e^{ikx_l}$: \begin{split} \frac{A_n-A_{n-1}}{\Delta t} + c \frac{A_ne^{ik\Delta x} - A_ne^{-ik\Delta x}}{2\Delta x} = 0 \end{split}
\begin{split}\downarrow\end{split}
\begin{split} A_n = A_{n-1} - \frac{c\Delta t}{2\Delta x} A_n (e^{ik\Delta x} - e^{-ik\Delta x}) \end{split}
\begin{split}\downarrow\end{split}
\begin{split} 1 = \frac{A_{n-1}}{A_n} + \frac{c\Delta t}{\Delta x}i\sin(k\Delta x) \end{split}
With some rearranging: \begin{split} \frac{A_n}{A_{n-1}} = \frac{1}{1 - \frac{c\Delta t}{\Delta x}i\sin(k\Delta x)} \end{split}
The denominator of this expression has a minimum magnitude of 1 and is otherwise greater than 1, so the amplification factor is always less than 1. This indicates that the numerical method is unconditionally stable.
Modified Equation Analysis¶
In addition to the expansions shown in the FTCS result, the following expansion is used to go backwards in time: \begin{split} u_l^{n-1} &= u - \Delta t u_t + \frac{\Delta t^2}{2}u_{tt} + \cdots \\ \end{split}
The discretized PDE is then: \begin{split} \frac{\Delta t u_t - \frac{\Delta t^2}{2}u_{tt}}{\Delta t} + c \frac{2\Delta x u_x + \frac{\Delta x^3}{3}u_{xxx}}{2\Delta x}=0 \end{split}
Again, $u_{tt} = c^2u_{xx}$.
\begin{split} u_t + cu_x = \frac{\Delta t}{2}u_{tt} - \frac{c\Delta x^2}{6}u_{xxx} \\ \end{split}
\begin{split}\downarrow\end{split}
\begin{split} u_t + cu_x = \frac{c^2 \Delta t}{2}u_{xx} - \frac{c\Delta x^2}{6}u_{xxx} \end{split}
The numerical diffusion coefficient is $\frac{c^2 \Delta t}{2}$, which is positive and thus diffusive unlike for the FTCS scheme. The numerical dispersion coefficient is $\frac{c\Delta x^2}{6}$.
This unconditional stability can be demonstrated numerically. Note that the implicit timestep requires solution of a linear system. The method of successive over-relaxation (SOR) will be used to solve this system.
In the solution, the diffusive effect of the numerical method can be observed - despite the original PDE having no diffusion or viscous term, the solution shows smearing.
def solverSOR(A, b, relax = 1.0, maxitr: int = 10000, tol: float = 1e-7):
M, N = A.shape
x = np.zeros_like(b, dtype=float)
for itr in range(maxitr):
x_prev = x.copy()
for j in range(N):
new_value = (b[j] - np.dot(A[j, :j], x[:j]) - np.dot(A[j, j+1:], x[j+1:])) / A[j, j]
x[j] = new_value * relax + (1 - relax) * x[j]
if np.max(np.abs(x - x_prev)) < tol:
# print(f"Converged in {itr + 1} iterations")
break
else:
print(f"Warning: did not converge within {maxitr} iterations")
return x
solution_BTCS = np.zeros((NX, NT))
solution_BTCS[:, 0] = u0
A = np.zeros((NX, NX))
b = np.zeros(NX)
for tdx in range(1, NT):
# build coefficient matrix
for xdx in range(NX):
lp1 = (xdx + 1) % NX
lm1 = (xdx - 1) % NX
A[xdx, lm1] = -C / (2 * dx)
A[xdx, xdx] = 1 / dt
A[xdx, lp1] = C / (2 * dx)
b[xdx] = solution_BTCS[xdx, tdx-1] / dt
solution_BTCS[:, tdx] = solverSOR(A, b, relax=0.5)
plt.figure(figsize=(7, 7))
plt.subplot(2, 1, 1)
for tdx in [0, 50, 100, 150]:
plt.plot(xarr, solution_BTCS[:, tdx], label=f"{Tarr[tdx]:.2f} s")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.title("Implicit Timestep, Centered Spatial Difference")
plt.subplot(2, 1, 2)
plt.imshow(solution_BTCS.T, origin="lower", aspect="auto", extent=(0, L, 0, T))
plt.colorbar(label="u(x)")
plt.xlabel("x")
plt.ylabel("t")
Text(0, 0.5, 't')
3: First-order Upwinding¶
Assuming $c>0$, then the discretization is: \begin{split} \frac{u_l^{n+1}-u_l^n}{\Delta t} + c \frac{u_{l}^n-u_{l-1}^n}{\Delta x} = 0 \end{split}
Stability Analysis¶
Assuming a solution of form $u_l^n = A_n e^{ikx_l}$: \begin{split} \frac{A_{n+1}e^{ik\Delta x} - A_ne^{ik\Delta x}}{\Delta t} + c \frac{A_ne^{ikx_l} - A_ne^{ik(x_l-\Delta x)}}{\Delta x} = 0 \end{split}
Divide out the exponential term and group terms:
\begin{split} A_{n+1} = A_n - \frac{c\Delta t}{\Delta x} A_n (1 - e^{-ik\Delta x}) = 0 \end{split}
Using $e^{-ik\Delta x} = \cos(k\Delta x) - i\sin(k\Delta x)$: \begin{split} \frac{A_{n+1}}{A_n} = 1 - \frac{c\Delta t}{\Delta x} + \frac{c\Delta t}{\Delta x}\cos(k\Delta x) - \frac{c\Delta t}{\Delta x}i\sin(k\Delta x) \end{split}
The amplification ratio can then be found as follows. For brevity, let $\frac{c\Delta t}{\Delta x} \equiv CFL$.
\begin{split} |G|^2 &= (1 - CFL + CFL \cos(k\Delta x))^2 + CFL^2 \sin^2(k\Delta x) \\ &= (1-CFL)^2 + 2(1-CFL)CFL\cos(k\Delta x) + CFL^2 \cos^2(k\Delta x) + CFL^2 \sin^2(k\Delta x) \\ &= 1 - 2CFL + CFL^2 + 2(CFL - CFL^2)\cos(k\Delta x) + CFL^2 \\ &= 1 - 2CFL + 2CFL^2 + 2(CFL - CFL^2)\cos(k\Delta x) \\ &= 1 - 2CFL(1-CFL)(1-\cos(k\Delta x)) \end{split}
Note that $(1-\cos(k\Delta x)) \geq 0 \space \forall \space x$. So, for the amplification ratio to be less than 1, we need $(1-CFL) \geq 0$. This gives us the stability criterion of $CFL \leq 1$, or $\frac{c \Delta t}{\Delta x} \leq 1$.
Modified Equation Analysis¶
The same expansion terms as previously described are used again.
Plugging into the discretized upwind PDE, the modified equation is then: \begin{split} \frac{\Delta t u_t + \frac{\Delta t^2}{2}u_{tt}}{\Delta t} + c \frac{\Delta x u_x - \frac{\Delta x^2}{2}u_{xx} + \frac{\Delta x^3}{6}u_{xxx}}{\Delta x} = 0 \end{split}
\begin{split}\downarrow\end{split}
\begin{split} u_t + \frac{\Delta t}{2}u_{tt} + c \big(u_x - \frac{\Delta x}{2}u_{xx} + \frac{\Delta x^2}{6}u_{xxx}\big) = 0 \end{split}
Collecting terms, and using $u_{tt}=c^2u_{xx}$, we get: \begin{split} u_t + cu_x = \bigg(\frac{c\Delta x}{2} - \frac{c^2 \Delta t}{2}\bigg) u_{xx} - \frac{c \Delta x^2}{6}u_{xxx} \end{split}
The numerical diffusion coefficient is $\big(\frac{c\Delta x}{2} - \frac{c^2 \Delta t}{2}\big)$ and the numerical dispersion coefficient is $-\frac{c \Delta x^2}{6}$.
solution_upwind = np.zeros((NX, NT))
solution_upwind[:, 0] = u0
for tdx in range(0, NT-1):
for xdx in range(NX):
lm1 = (xdx - 1) % NX
solution_upwind[xdx, tdx+1] = (
solution_upwind[xdx, tdx]
+ dt * (-C * (solution_upwind[xdx, tdx] - solution_upwind[lm1, tdx])/dx)
)
plt.figure(figsize=(7, 7))
plt.subplot(2, 1, 1)
for tdx in [0, 50, 100, 150]:
plt.plot(xarr, solution_upwind[:, tdx], label=f"{Tarr[tdx]:.2f} s")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.title("First-order Upwinding")
plt.subplot(2, 1, 2)
plt.imshow(solution_upwind.T, origin="lower", aspect="auto", extent=(0, L, 0, T))
plt.colorbar(label="u(x)")
plt.xlabel("x")
plt.ylabel("t")
Text(0, 0.5, 't')
4: Lax-Friedrichs¶
The discretized PDE is: \begin{split} \frac{u_l^{n+1} - \frac{1}{2}(u_{l+1}^n + u_{l+1}^n)}{\Delta t} + c \frac{u_{l+1}^n - u_{l+1}^n}{2\Delta x} = 0 \end{split}
Stability Analysis¶
Assuming a solution of form $u_l^n = A_n e^{ikx_l}$: \begin{split} \frac{A_{n+1}e^{ik\Delta x} - \frac{1}{2}(A_ne^{ik(x+\Delta x)} + A_ne^{ik(x-\Delta x)})}{\Delta t} + c \frac{A_ne^{ik(x+\Delta x)} - A_ne^{ik(x-\Delta x)}}{2\Delta x} = 0 \end{split}
\begin{split}\downarrow\end{split}
\begin{split} A_{n+1} = \frac{1}{2}\big( A_n e^{ik\Delta x} + A_n e^{-ik\Delta x} \big) + \frac{c\Delta t}{2\Delta x} \big( A_ne^{ik\Delta x} - A_ne^{-ik\Delta x} \big) = 0 \end{split}
\begin{split}\downarrow\end{split}
\begin{split} \frac{A_{n+1}}{A_n} = \frac{1}{2}\big( e^{ik\Delta x} + e^{-ik\Delta x} \big) + \frac{c\Delta t}{2\Delta x} \big( e^{ik\Delta x} - e^{-ik\Delta x} \big) = 0 \end{split}
Using Euler's identity ($e^{ik\Delta x} = \cos(k\Delta x) + i \sin(k\Delta x)$): \begin{split} \frac{A_{n+1}}{A_n} = \cos(k\Delta x) - \frac{c\Delta t}{\Delta x}i\sin(k\Delta x) \end{split}
The amplification ratio is then: \begin{split} \bigg|\frac{A_{n+1}}{A_n}\bigg|^2 &= \cos^2(k\Delta x) + \bigg(\frac{c\Delta t}{\Delta x}\bigg)^2\sin^2(k\Delta x) \\ &= 1 + \bigg(\big(\frac{c\Delta t}{\Delta x}\big)^2 - 1\bigg)\sin^2(k\Delta x) \end{split}
For stability, we need $\big((\frac{c\Delta t}{\Delta x})^2 - 1\big)\sin^2(k\Delta x) \leq 0$. The $\sin^2$ term is always positive, so that leaves us with $\frac{c\Delta t}{\Delta x} \leq 1$ as the stability criterion.
Modified Equation Analysis¶
The same expansion terms as previously described are used again.
Collecting some terms first: \begin{split} u_{l+1}^n + u_{l-1}^n &= 2u + \Delta x^2 u_{xx} \\ u_{l+1}^n - u_{l-1}^n &= 2\Delta x u_x + \frac{\Delta x^3}{3}u_{xxx} \\ \end{split}
Plugging into the Lax-Friedrichs discretized PDE: \begin{split} \frac{\Delta t u_t + \frac{\Delta t^2}{2}u_{tt} - \frac{\Delta x^2}{2}u_{xx}}{\Delta t} + c \frac{2\Delta x u_x + \frac{\Delta x^3}{3}u_{xxx}}{2\Delta x} = 0 \end{split}
\begin{split}\downarrow\end{split}
\begin{split} u_t + \frac{\Delta t}{2}u_{tt} - \frac{\Delta x^2}{2\Delta t}u_{xx} + c\big(u_x + \frac{\Delta x^2}{6}u_{xxx}\big)=0 \end{split}
\begin{split}\downarrow\end{split}
\begin{split} u_t + cu_x = -\frac{\Delta t}{2}u_{tt} + \frac{\Delta x^2}{2\Delta t}u_{xx} - c\frac{\Delta x^2}{6}u_{xxx} \end{split}
The numerical diffusion coefficient is then found to be $\frac{\Delta x^2}{2\Delta t} - \frac{c^2 \Delta t}{2}$, and the nuermical dispersion is found to be $- c\frac{\Delta x^2}{6}u_{xxx}$.
solution_lf = np.zeros((NX, NT))
solution_lf[:, 0] = u0
for tdx in range(0, NT-1):
for xdx in range(NX):
lp1 = (xdx + 1) % NX
lm1 = (xdx - 1) % NX
solution_lf[xdx, tdx+1] = (
(solution_lf[lp1, tdx] + solution_lf[lm1, tdx]) / 2
+ dt * (-C * (solution_lf[lp1, tdx] - solution_lf[lm1, tdx])/(2*dx))
)
plt.figure(figsize=(7, 7))
plt.subplot(2, 1, 1)
for tdx in [0, 50, 100, 150]:
plt.plot(xarr, solution_lf[:, tdx], label=f"{Tarr[tdx]:.2f} s")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.title("Lax-Friedrichs")
plt.subplot(2, 1, 2)
plt.imshow(solution_lf.T, origin="lower", aspect="auto", extent=(0, L, 0, T))
plt.colorbar(label="u(x)")
plt.xlabel("x")
plt.ylabel("t")
Text(0, 0.5, 't')
5: Lax-Wendroff¶
Here, the Lax-Wendroff method is applied to the first-order upwind scheme, for which we know the numerical diffusion is $\frac{c\Delta x}{2}(1-CFL)$. This numerical diffusion is "removed" from the discretized equation as an anti-diffusive term. That looks like: \begin{split} u_t + cu_x = -\frac{c\Delta x}{2}(1-CFL)u_{xx} \end{split} which is then discretized as follows: \begin{split} u_l^{n+1} = u_l^n - \frac{CFL}{2}(u_{l}^n - u_{l-1}^n) - \frac{CFL(1-CFL)}{2}(u_{l+1}^n - 2u_l^n + u_{l-1}^n) \end{split}
Stability Analysis¶
Assuming $u_l^n=A_ne^{ikx_l}$: \begin{split} A_{n+1}e^{ikx_l} = A_ne^{ikx_l} - \frac{CFL}{2}(A_ne^{ikx_l} - A_ne^{ik(x_l-\Delta x)}) - \frac{CFL(1-CFL)}{2}A_n(e^{ik(x_l+\Delta x)}-2e^{ikx_l}+e^{ik(x_l-\Delta x)}) \end{split}
This simplifies to: \begin{split} G = \frac{A_{n+1}}{A_n} = 1 - CFL(1-e^{-ik\Delta x}) - \frac{CFL(1-CFL)}{2}(2\cos(k\Delta x)-2) \end{split}
The real part of the above expression is $\text{Re}(G) = 1 + CFL^2(\cos(k\Delta x)-1)$, and the imaginary is $\text{Im}(G) = -CFL\sin(k\Delta x)$.
Therefore: $G = 1 + CFL^2(\cos(k\Delta x)-1) -iCFL\sin(k\Delta x)$.
Expanding this expression to obtain the magnitude $|G|$ for stability analysis results in disgusting math. The easier way is to consider limiting wavenumber cases to simplify the expression first before expanding.
Starting with $k=\pi/2$: \begin{split} \big|G\big|_{k=\pi/2} = (1-CFL)^2 + CFL^2 \leq 1 \\ 1 - CFL^2 + CFL^4 \leq 1 \\ CFL^4 - CFL^2 \leq 0 \\ CFL^2 - 1 \leq 0 \end{split}
This can also be done with $k=\pi$. Both will result in an expression that finds the stability limit to be $CFL \leq 1.0$.
Modified Equation Analysis¶
By definition, the Lax-Wendroff scheme on the first-order upwinding method removes the numerical diffusion.
The only term left is the numerical dispersion term, which is $-\frac{c \Delta x^2}{6}$.
This is the same value as for the pure first-order upwinding scheme.
In the numerical simulation below, the effect of this dispersion can be seen as the solution starts to "ring" around sharp transitions.
It should be noted that the cancellation of the numerical diffusion by the Lax-Wendroff scheme introduces higher order errors. In this case, additional third order error is introduced, but this does not affect the fact that the scheme is first order accurate in time and second order in space.
solution_lw = np.zeros((NX, NT))
solution_lw[:, 0] = u0
for tdx in range(0, NT-1):
for xdx in range(NX):
lp1 = (xdx + 1) % NX
lm1 = (xdx - 1) % NX
solution_lw[xdx, tdx+1] = (
solution_lw[xdx, tdx]
+ dt * (-C * (solution_lw[xdx, tdx] - solution_lw[lm1, tdx])/dx)
- CFL * (1 - CFL) / 2 * (solution_lw[lp1, tdx] - 2 * solution_lw[xdx, tdx] + solution_lw[lm1, tdx])
)
plt.figure(figsize=(7, 7))
plt.subplot(2, 1, 1)
for tdx in [0, 50, 100, 150]:
plt.plot(xarr, solution_lw[:, tdx], label=f"{Tarr[tdx]:.2f} s")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.title("Lax-Wendroff (based on First-order Upwinding)")
plt.subplot(2, 1, 2)
plt.imshow(solution_upwind.T, origin="lower", aspect="auto", extent=(0, L, 0, T))
plt.colorbar(label="u(x)")
plt.xlabel("x")
plt.ylabel("t")
Text(0, 0.5, 't')
6 Crank-Nicholson¶
Following the Crank-Nicholson scheme, $u_t + cu_x = 0$ is discretized as: \begin{split} \frac{u_l^{n+1} - u_l^n}{\Delta t} + c \frac{ u_{l+1}^{n+1} - u_{l-1}^{n+1} + u_{l+1}^{n} - u_{l-1}^{n} }{4\Delta x} = 0 \end{split}
This can be rearranged to be: \begin{split} u_l^{n+1} + \frac{CFL}{4}\big(u_{l+1}^{n+1} - u_{l-1}^{n+1}\big) = u_l^n - \frac{CFL}{4}\big(u_{l+1}^{n} - u_{l-1}^{n}\big) \end{split}
Stability Analysis¶
Using $u_l^n=A_ne^{ikx_l}$: \begin{split} A_{n+1}e^{ikx_l} + \frac{CFL}{4}\big( A_{n+1}e^{ik(x_l+\Delta x)} - A_{n+1}e^{ik(x_l-\Delta x)} \big) = A_ne^{ikx_l} - \frac{CFL}{4}\big( A_{n}e^{ik(x_l+\Delta x)} - A_{n}e^{ik(x_l-\Delta x)} \big) \end{split}
Divide out $e^{ikx_l}$: \begin{split} A_{n+1} + \frac{CFL}{4}\big( A_{n+1}e^{ik\Delta x} - A_{n+1}e^{-ik\Delta x} \big) = A_n - \frac{CFL}{4}\big( A_{n}e^{ik\Delta x} - A_{n}e^{-ik\Delta x} \big) \end{split}
Using $e^{ik\Delta x} - e^{-ik\Delta x} = 2i\sin(k\Delta x)$, the following expression for amplification ratio can be obtained: \begin{split} \frac{A_{n+1}}{A_n} = \frac{1 - i\frac{CFL}{2}\sin(k\Delta x)}{1 + i\frac{CFL}{2}\sin(k\Delta x)} \end{split}
From this, it can be seen that $\big|\frac{A_{n+1}}{A_n}\big|=1$ always, meaning that the Crank-Nicholson method is unconditionally stable (but right on the edge). It can also be inferred that theoretically there should be no numerical diffusion due to the amplification ratio being 1 exactly.
Modified Equation Analysis¶
The same expansions as used previously are still used here. However, the $n+1$ terms must be handled using a 2D Taylor expansion, which is as follows: \begin{split} f(a+h,b+k) = f + hf_a + kf_b + \frac{h^2}{2}f_{aa} + \frac{k^2}{2}f_{bb} + hkf_{ab} + \cdots \end{split} Therefore: \begin{split} u_{l\pm 1}^{n+1} = u \pm \Delta x u_x + \Delta t u_t \pm \Delta x \Delta t u_{xt} + \frac{\Delta x^2}{2}u_{xx} + \frac{\Delta t^2}{2}u_{tt} \pm \frac{\Delta x^3}{6}u_{xxx} + \cdots \end{split}
Starting from the Crank-Nicholson discretization: \begin{split} u_{l+1}^{n+1} - u_{l-1}^{n+1} &= 2\Delta x u_x + \frac{\Delta x^3}{3} u_{xxx} + 2\Delta x \Delta t u_{xt} \\ u_{l+1}^{n} - u_{l-1}^{n} &= 2\Delta x u_x + \frac{\Delta x^3}{3} u_{xxx} \end{split} So the numerator of the fraction is: \begin{split} u_{l+1}^{n+1} - u_{l-1}^{n+1} + u_{l+1}^{n} - u_{l-1}^{n} = 4\Delta x u_x + \frac{2 \Delta x^3}{3}u_{xxx} + 2\Delta x \Delta t u_{xt} \end{split}
Substituting back into the discretized PDE, the modified equation is: \begin{split} \frac{\Delta t u_t + \frac{\Delta t^2}{2}u_{tt}}{\Delta t} + \frac{c}{4\Delta x} \bigg( 4\Delta x u_x + \frac{2 \Delta x^3}{3}u_{xxx} + 2\Delta x \Delta t u_{xt} \bigg) = 0 \\ u_t + \frac{\Delta t}{2}u_{tt} + c \bigg( u_x + \frac{\Delta^2}{6} u_{xxx} + \frac{\Delta t}{2}u_{xt} \bigg) = 0 \end{split}
Since $u_t=-cu_x$, $u_{tt} = -cu_{xt}$. As such: \begin{split} u_t + cu_x = -\frac{c\Delta x^2}{6}u_{xxx} \end{split}
Therefore the numerical dispersion coefficient is $-\frac{c\Delta x^2}{6}$, and there is no numerical diffusion.
solution_cn = np.zeros((NX, NT))
solution_cn[:, 0] = u0
A = np.zeros((NX, NX))
b = np.zeros(NX)
for tdx in range(1, NT):
# build coefficient matrix
for xdx in range(NX):
lp1 = (xdx + 1) % NX
lm1 = (xdx - 1) % NX
A[xdx, lm1] = -CFL/4
A[xdx, xdx] = 1
A[xdx, lp1] = CFL/4
b[xdx] = solution_cn[xdx, tdx-1] - (CFL / 4) * (solution_cn[lp1, tdx-1] - solution_cn[lm1, tdx-1])
solution_cn[:, tdx] = solverSOR(A, b, relax=0.8)
plt.figure(figsize=(7, 7))
plt.subplot(2, 1, 1)
for tdx in [0, 50, 100, 150]:
plt.plot(xarr, solution_cn[:, tdx], label=f"{Tarr[tdx]:.2f} s")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.title("Crank-Nicholson")
plt.subplot(2, 1, 2)
plt.imshow(solution_cn.T, origin="lower", aspect="auto", extent=(0, L, 0, T))
plt.colorbar(label="u(x)")
plt.xlabel("x")
plt.ylabel("t")
Text(0, 0.5, 't')
7 Spectral Method¶
The Fourier transform of a function $f(x)$ is: \begin{split} \mathcal{F}(f(x)) = \hat{f}(k) = \int f(x) e^{-i k x} dx \end{split}
The inverse Fourier transform is: \begin{split} \mathcal{F}^{-1}(\hat{f}(k)) = f(x) = \int \hat{f}(k) e^{i k x} dk \end{split}
Starting with the inverse Fourier transform, the derivative in Fourier space is: \begin{split} \frac{d}{dx}f(x) = \int \hat{f}(k) \frac{d}{dx}e^{ikx} dk = \int ik \hat{f}(k) e^{ikx} dk \end{split}
Taking the Fourier transform of the above, we find that $\frac{d}{dx}\hat{f} = ik \hat{f}(k)$.
For convienience, the FFT module in NumPy will be used to perform the Fourier transform and obtain the wavenumbers.
An RK4 integrator is used to help with stability.
Note that even though mathematically the FFT spectral method is exact, in practice due to truncated Fourier series there is frequency truncation. The loss of the high frequency information results in "ringing" (Gibbs phenomenon) near sharp corners. It may look like dispersion but is not caused by uneven wavespeeds.
solution_spectral = np.zeros((NX, NT))
solution_spectral[:, 0] = u0
k = np.fft.fftfreq(NX, d=dx/(2*np.pi))
fn_u_x = lambda x: np.fft.ifft(1j * k * np.fft.fft(x)).real
fn_u_t = lambda x: -C * fn_u_x(x)
for tdx in range(0, NT-1):
k1 = fn_u_t(solution_spectral[:, tdx])
k2 = fn_u_t(solution_spectral[:, tdx] + k1 * dt / 2)
k3 = fn_u_t(solution_spectral[:, tdx] + k2 * dt / 2)
k4 = fn_u_t(solution_spectral[:, tdx] + k3 * dt)
solution_spectral[:, tdx + 1] = solution_spectral[:, tdx] + (k1 + 2*k2 + 2*k3 + k4) * dt / 6
plt.figure(figsize=(7, 7))
plt.subplot(2, 1, 1)
for tdx in [0, 50, 100, 150]:
plt.plot(xarr, solution_spectral[:, tdx], label=f"{Tarr[tdx]:.2f} s")
plt.legend()
plt.xlabel("x")
plt.ylabel("u(x)")
plt.title("Spectral Method")
plt.subplot(2, 1, 2)
plt.imshow(solution_spectral.T, origin="lower", aspect="auto", extent=(0, L, 0, T))
plt.colorbar(label="u(x)")
plt.xlabel("x")
plt.ylabel("t")
Text(0, 0.5, 't')
