Higher Order PDE using Finite Difference

Question

How to approximate higher-order partial differential equation using finite difference method?

$$\frac{\partial^{2} y}{\partial t^{2}}+\frac{\partial^{4} y}{\partial x^{4}}=0$$

Answer 1

This is an equation of a dynamic Euler–Bernoulli beam. It can be simply approximated explicitly using standard second and fourth derivative approximation. $$ \frac{\partial^2 u}{\partial t^2} \sim \frac{u^{n+1} - 2u^{n} + u^{n-1}}{\Delta t^2}\\ \frac{\partial^4 u}{\partial x^4} \sim \frac{u_{m-2} - 4u_{m-1} + 6u_m -4u_{m+1}+u_{m+2}}{\Delta x^4} $$ Both formulas have truncation error of second order $O(\Delta t^2)$ and $O(\Delta x^2)$ due to symmetry. There are plenty of methods to derive the approximation for the fourth derivative. The simplest one would be applying second derivative formula twice $$ \frac{\partial^4 u}{\partial x^4} = \frac{\partial^2}{\partial x^2}\frac{\partial^2 u}{\partial x^2} \sim \frac{1}{\Delta x^2}\left[ \left(\frac{u_{m+2} - 2u_{m+1}+u_{m}}{\Delta x^2}\right) -2\left(\frac{u_{m+1} - 2u_m+u_{m+1}}{\Delta x^2}\right) +\left(\frac{u_{m} - 2u_{m-1}+u_{m-2}}{\Delta x^2}\right) \right] $$ The coefficient are the same as in binomial theorem for $(a-b)^4$, since $$ \frac{u_{m-2} - 4u_{m-1} + 6u_m -4u_{m+1}+u_{m+2}}{\Delta x^4} = \left(\frac{T_{+1/2} - T_{-1/2}}{\Delta x}\right)^4 u_m $$ where $T_k$ are shift operators: $T_k u_m \equiv e^{k\Delta x \frac{\partial}{\partial x}}u_m = u_{m+k}$.

So now we can construct a $O(\Delta t^2, \Delta x^2)$ explicit scheme of form $$ \frac{u^{n+1}_m - 2u^n_m + u^n_{m+1}}{\Delta t^2} + K\frac{u^n_{m-2} - 4u^n_{m-1} + 6u^n_m -4u^n_{m+1}+u^n_{m+2}}{\Delta x^4} = 0. $$ Here $K = \frac{EI}{\rho} = 1$ for your equation. Let's investigate whether it is stable using von Neumann method. Since the equation is a difference equation with constant coefficients we can seek a solution in form $u^n_m = \lambda^n e^{i\alpha m}$. Plugging that into the equation and canceling $\lambda^n e^{i\alpha m}$ gives $$ \lambda -2 + \frac{1}{\lambda} + \sigma \left(e^{i\alpha/2} - e^{-i\alpha/2}\right)^4 = 0 $$ where $\sigma$ is the Courant's number for your problem, specifically $\sigma = \frac{K \Delta t^2}{\Delta x^4}$. Changing $e^{i\alpha/2} - e^{-i\alpha/2} = 2\sin \frac{\alpha}{2}$ gives a quadratic equation with real coefficients for $\lambda$ $$ \lambda^2 -2\lambda + 1 + 16\sigma\lambda \sin^4\frac{\alpha}{2} = 0. \tag{*} $$ For a scheme to be stable every $\lambda$ for any $\alpha$ should not exceed 1 by its absolute value. If the $(*)$ equation has two real roots they could not both satisfy $|\lambda_{1,2}| < 1$ since their product is $\lambda_1 \lambda_2 = 1$ due to Vieta's formulas. So there left two possibilities. Either $\lambda_{1,2}$ is a complex-conjugated pair or the roots are real and $\lambda_1 = \lambda_2 = \pm 1$. So the discriminant of the quadratic should be negative or equal to zero. $$ \frac{D}{4} = \left(1 - 8\sigma \sin^2 \frac{\alpha}{2}\right)^2 - 1 = -16\left(\sigma \sin^2 \frac{\alpha}{2}\right)\left(1 - 4\sigma \sin^2 \frac{\alpha}{2}\right) \leqslant 0. $$ So for the scheme to be stable the following condition must hold for every $\alpha$ $$ 4\sigma \sin^2 \frac{\alpha}{2} \leqslant 1 $$ which would hold iff $$ \frac{K\Delta t^2}{\Delta x^4} \equiv \sigma \leqslant \frac{1}{4}\\ \Delta t \leqslant \frac{\sqrt{K}h^2}{2}. $$