I'm trying to solve a nonlinear partial differential equation \begin{equation} L(u_{xxtt},u_{xx}u_{tt},u_{xt}^2,u_{xt},u_{tt})=0 \end{equation} using finite difference methods. In order to remove the second derivatives with respect to time, I substitute \begin{equation} v := u_t \end{equation} such that it becomes \begin{equation} L(v_{xxt},u_{xx}v_{t},v_{x}^2,v_{x},v_{t}) = 0 \end{equation} So in my algorithm, I first implicitly calculate $v$ at the next time step using central difference approximations (in order to ensure that local truncation error is second order both in time and space), i.e. \begin{align} (v_i{}^m)_{xx} & = (v^m_{i+1} - 2 v_i{}^m + v^m_{i-1})/\Delta x^2 \\ (v_i{}^m)_x & = (v^m_{i+1}-v^m_{i-1}) / 2\Delta x \\ (v_i{}^m)_t & = (v_i^{m+1} - v_i{}^{m-1})/2 \Delta t \end{align} Once I've calculated $v_i^{m+1}$, I calculate $u_i^{m+1}$ as follows \begin{equation} v_i{}^m = (u_i{}^m)_t = \frac{u^{m+1}_i - u_i^{m-1}}{2 \Delta t} \implies u^{m+1} = 2 \Delta t v_i{}^m + u_i^{m-1} \tag{1} \end{equation} I would like to determine the stability of this program, and any tips on how to do that are greatly appreciated.
I first wanted to determine the stability using von Neumann's method, but I don't think this is possible (but please correct me if I'm wrong). Let me explain:
Based on the answer given to my previous question, I've decided to "freeze" $u_{xx}$ and $v_x$ in order to linearize the PDE \begin{equation} L(v_{xxt},v_{x},v_{t}) = 0 \end{equation} However, if I substitute equation $(1)$ into my finite difference operator, then I will get terms involving $u^{m+2}, u^{m+1}, u^m, u^{m-1}, u^{m-2}$. Thus, if I now substitute in the following (as per usual for the von Neumann stability analysis) \begin{equation} u_i{}^m = A^m \exp(\hat{i} k i \Delta x) \; , \; \hat{i}:= \sqrt{-1} \end{equation} then I will end up with a fifth order polynomial in $A$ \begin{equation} aA^4 + bA^3 + cA^2 + dA + e = 0 \end{equation} This means it will not be possible to solve for $A$ and thus I cannot determine the amplification factor. Hence, unless I've missed something, it's not possible to use the von Neumann analysis to determine stability of the scheme.
