Stability analysis of coupled ordinary differential equations

Question

Given a forward-in-time approximation I have the coupled equations: $$ \frac{T^{(n+1)} - T^{(n)}}{\Delta t} = x T^{(n)} - y h^{(n)} \\ \frac{h^{(n+1)} - h^{(n)}}{\Delta t} = -z h^{(n)} - \alpha T^{(n)} $$ where $x, y, z$ and $\alpha$ are constants, I see from my simulation that the solution is damping. But how can I use Von Neumann analysis to find the amplication factor $A$?

Answer 1

You can rewrite your finite difference method into the form

\begin{align}\left[\begin{matrix} T^{t+1} \\ h^{t+1}\end{matrix} \right]=\left[\begin{matrix}A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix} \right] \left[\begin{matrix} T^n\\h^n\end{matrix} \right] = \left[\begin{matrix}A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix} \right]^n\left[\begin{matrix} T^0\\h^0\end{matrix} \right] \end{align}

Rewritten this way, you will find that the matrix $A$ becomes $$A=I+\Delta tJ,$$ where $I$ is the identity matrix and $J$ is the jacobian of the RHS of your original system. $A$ can be thought of as the amplification matrix, whose eigenvalues determine the stability of the numerical system. To ensure absolute stability, you need to find a region in the complex vector space such that the spectral radius $\rho(A)<1$ (i.e. the dominant eigenvalue is less than one).

Of course, since this is a system of ODE's, you need to be extra careful. Stability and consistency may not be enough to ensure an accurate solution. You'll also need to check for Stiffness.