In the Euler-Bernoulli model the equations of motion of a vibrating beam are given by $$ \mu X_{tt}+EIX_{zzzz}=0, $$ where $X(z, t)$ is the deflexion at height $z$ and at time $t$. We discretize this function in the following way $$ X(z, t) \rightarrow X(z_j, t_n)\rightarrow X(z_0+jh, n\delta t)\rightarrow X_{j,n}, $$ where $h$ is the grid size and $\delta t$ is the integration time-step. The time derivative is discretized by the following approximation $$ X_{tt}(z, t)\approx \frac{X_{j,n+1}+X_{j, n-1}-2X_{j,n}}{\delta t^2}, $$ which has an error of order $\mathcal{O}(\delta t^2)$. The central finite-difference discretization of $X_{zzzz}$ used at interior points of the grid is given by $$ X_{zzzz}(z, t)\approx\dfrac{X_{j+2, n}-4X_{j+1, n}+6X_{j, n}-4X_{j-1, n}+X_{j-2, n}}{h^4}, $$ whose error is of the order of $\mathcal{O}(h^2)$. A simple and straight forward explicit method for the numerical time integration of interior points of the grid is given by $$ X_{j, n+1}= 2X_{j, n}-X_{j, n-1}+\dfrac{\delta t^2EI}{h^4\mu}(X_{j+2, n} -4X_{j+1, n}+6X_{j, n}-4X_{j-1, n}+X_{j-2, n}) $$ We use clamped boundary conditions at $z=0$, $X(0, t)=X_z(0, t)=0$ and free boundary conditions at $z=\ell_0$, $X_{zz}(\ell_0, t)=X_{zzz}(\ell_0, t)=0$.
We next analyse the stability of this integration scheme using the von Neumann stability criterion. Suppose $X_{j,n}=g^ne^{ikz_j}$, where $z_j=z_0+jh$. We obtain then \begin{align} g&=2-1/g+\dfrac{\delta t^2EI}{h^4\mu}\left[e^{i2kh}-4e^{ikh}+6-4e^{-ikh}+e^{-2ikh}\right]\\ &=2-1/g+\dfrac{2\delta t^2EI}{h^4\mu} \left[3-4\cos(kh)+\cos(2kh)\right] \end{align} with the notation $$ \beta = 1+\alpha\left[3-4\cos(kh)+\cos(2kh)\right], $$ where $\alpha=\dfrac{\delta t^2EI}{h^4\mu}$. We then find that the equation for $g$ is given by $$ g^2-2\beta g+1=0. $$ The solutions are given by $$ g_\pm= \beta\pm\sqrt{\beta^2-1}. $$ Hence, if $\beta>1$, this method of integration is unstable. One can show that $\beta>1$. Therefore, we have to find another integration method. We now modify the previous explicit method and implement an implicit method for the numerical time integration. It is given by \begin{align} X_{i, n+1} &= 2X_{i, n}-X_{i, n-1}+\dfrac{\delta t^2EI}{2h^4\mu} \left[X_{i+2, n}+X_{i+2, n+1} -4(X_{i+1, n}+X_{i+1, n+1})\right.\\ &\left.+6(X_{i, n}+X_{i, n+1})-4(X_{i-1, n}+X_{i-1, n+1})+X_{i-2, n}+X_{i-2, n+1}\right] \end{align} With the substitution $X_{j,n}=g^ne^{ikz_j}$, we find $$ \kappa g-\gamma g+1=0, $$ where \begin{align} \kappa &=1-\alpha[3-4\cosh(kh)+\cosh(2kh)],\\ \gamma &=2+\alpha[3-4\cosh(kh)+\cosh(2kh)]. \end{align} The solution for the quadratic equation in $g$ is $$ g_\pm=\frac{\gamma\pm\sqrt{\gamma^2-4\kappa}}{2\kappa}. $$ One can verify that $|g_\pm|>1$, hence this implies this method is also unstable. Could anyone provide a stable (finite difference or finite element) numerical integration method for the Euler-Bernoulli equations?