# Stability analysis of coupled ordinary differential equations

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$?

• can someone explain von neumann stability analysis for system of equation Dec 8 '16 at 14:27

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.

• good answer. This may also explain the damping via A-stability. However, this is different from von Neumann stability analysis.
– Jan
Jan 28 '14 at 20:42
• @jan: indeed this is not von neumann analysis, for the equation itself isn't a PDE, but an ODE.
– Paul
Jan 28 '14 at 20:53
• Ouch, sounds like I need to learn me some more numerical methods before I attempt this! We've not covered stiffness or eigenvalues. Jan 28 '14 at 21:38