Solving fourth order diffusion type equation using Finite Difference method

86 Views Asked by At

I am trying to solve a fourth order partial differential equation $$ {\partial u\over \partial t} = i{\partial^4 u\over \partial x^4} $$ where $i=\sqrt{-1}$ and periodic boundary conditions.

I tried using first order in time and 2nd order in space finite difference scheme, even second order in time and space also. But I am unable to formulate a FD scheme with is stable. I am trying to check stability using von-Neumann stability criteria. Could somebody explain how to proceed in this case?

1

There are 1 best solutions below

0
On

Let's try the second order accurate central finite difference for the spatial variable and Backward Euler for the temporal discretization. Then, your discretized PDE reads: $$\frac{u_j^{k+1} - u_j^k}{\Delta t} = i \frac{u_{j-2}^{k+1} - 4u_{j-1}^{k+1} + 6 u_{j}^{k+1} -4 u_{j+1}^{k+1} + u_{j+2}^{k+1} }{\Delta x^4}$$ Now make the usual Ansatz for a von-Neumann stability analysis $$u_j^k = e^{ikx_j}$$ and get your hands dirty: \begin{align} u_j^{k+1} &= e^{ikx_j} + \underbrace{r}_{i \Delta t / \Delta x^4} \Big( e^{i(k+1) (x_j - 2 \Delta x)} - 4 e^{i(k+1) (x_j - \Delta x)} + 6e^{i(k+1)x_j} \\ & \quad \quad \quad \quad \quad \quad \quad - 4 e^{i(k+1) (x_j + \Delta x)} + e^{i(k+1) (x_j + 2\Delta x)} \Big) \\ &= e^{ikx_j} + r \Big( e^{ik x_j } e^{i (x_j - 2 \Delta xk - 2 \Delta x)} - 4 e^{ikx_j} e^{i(x_j -k\Delta x - \Delta x) } + 6e^{ik x_j } e^{ix_j} \\ &\quad \quad \quad \quad \quad - 4 e^{ik x_j } e^{i (x_j + k\Delta x + \Delta x)} + e^{ik x_j } e^{i(x_j + 2\Delta x + 2 \Delta x k)} \Big) \\ &= e^{ik x_j } \Big[ 1 + r \Big( e^{i (x_j - 2 \Delta xk - 2 \Delta x)} - 4 e^{i(x_j -k\Delta x - \Delta x) } + 6 e^{ix_j} \\ & \quad \quad \quad \quad \quad \quad \quad - 4 e^{i (x_j + k\Delta x + \Delta x)} + e^{i(x_j + 2\Delta x + 2 \Delta x k)} \Big) \Big] \end{align}

By Wolfram Alpha, \begin{align}& e^{i (x_j - 2 \Delta xk - 2 \Delta x)} - 4 e^{i(x_j -k\Delta x - \Delta x) } + 6 e^{ix_j} - 4 e^{i (x_j + k\Delta x + \Delta x)} + e^{i(x_j + 2\Delta x + 2 \Delta x k)} \\ =& 16 e^{ix_j} \sin^4\Big( 0.5\Delta x [k+1] \Big)\end{align}

And using Euler's formula $e^{ix_j} = \cos(x_j) + i \sin(x_j)$

\begin{align} &r \Big( e^{i (x_j - 2 \Delta xk - 2 \Delta x)} - 4 e^{i(x_j -k\Delta x - \Delta x) } + 6 e^{ix_j} - 4 e^{i (x_j + k\Delta x + \Delta x)} + e^{i(x_j + 2\Delta x + 2 \Delta x k)} \Big) \\ =& \underbrace{16 \frac{\Delta t}{\Delta x^4} \sin^4\Big( 0.5\Delta x [k+1] \Big)}_{=: \alpha(k)} \big[ i \cos(x_j) - \sin(x_j) \big] \end{align}

Thus, the absolute value of the amplification $$ |G| = \Big \vert 1 + r \Big( e^{i (x_j - 2 \Delta xk - 2 \Delta x)} - 4 e^{i(x_j -k\Delta x - \Delta x) } + 6 e^{ix_j} - 4 e^{i (x_j + k\Delta x + \Delta x)} + e^{i(x_j + 2\Delta x + 2 \Delta x k)} \Big) \Big \vert$$ is given by \begin{align} |G|^2 &= \alpha(k)^2 \cos(x_j)^2 + (1 - \alpha(k) \sin(x_j)^2) \\ &= \alpha(k)^2 \cos(x_j)^2 + \alpha(k)^2 \sin(x_j)^2+ 1 - 2 \alpha(k) \sin(x_j) \\ &= \alpha(k)^2 + 1 - 2 \alpha(k) = (\alpha(k) - 1)^2 \end{align} So to have stability $|G| < 1$, we require \begin{align} \vert \alpha(k) - 1 \vert < 1 \end{align} Note that it holds $\forall k \in \mathbb{N}_0, \Delta t > 0, \Delta x >0$: $$0 < \alpha (k) \leq 16 \frac{\Delta t}{\Delta x^4} $$ Thus if you guarantee that $$ 16 \frac{\Delta t}{\Delta x^4} -1 < 1$$ you have stability.