# How to estimate stability and stiffness of a system of coupled ODEs?

I'm running into issues with Python/Julia ODE solvers requiring prohibitively small timesteps to evolve a system of 4 coupled ODEs (the order of magnitude of the state variables and time unit span ~40-50 orders of magnitude). Interestingly if I use a simple explicit Euler approach then I can get past the stiffness at early times (where the general adaptive solvers get stuck) and roughly match theoretical expectations at later times.

Apart from the technical coding details, I have noticed that it is common to talk about the mathematical stability and stiffness of ODE systems. Is there a standard way to calculate the expected stability and stiffness given a set of coupled ODEs, initial conditions, integration timescales, input parameters, and maybe a solver algorithm?

For what it's worth, my system of 4 coupled ODEs looks like

$$\frac{dM_1}{dt} = g(t) - c(t,M_1,E)M_1(t) + w(t)M_3(t) - e(M_1,E)h(M_1,E)E(t)$$

$$\frac{dM_2}{dt} = c(t,M_1,E)M_1(t) - s(t)M_2(t) - w(t)M_3(t)$$

$$\frac{dM_3}{dt} = s(t)M_2(t)$$

$$\frac{dE}{dt} = g(t)v(t) - c(t,M_1,E)E(t) + b(t)w(t)M_3(t) - h(M_1,E)E(t)$$

where the various lettered functions are input parameters that depend on time and/or state variables. I know this isn't enough for doing stability/stiffness analysis but it's a start. I'd appreciate any pedagogical references, even for much simpler illustrative systems -- I vaguely remember this involves linear algebra and eigendecomposition but it's been so long since undergrad.

• Have you computed the eigenvalues of the Jacobian of the right hand side? Mar 30 at 15:13
• No but I would love to explore that idea further for different assumptions about the input parameters. Does that need to be done manually or is there a standard package in scipy that can return the Jacobian given a set of ODEs (or matrix of the RHS)? Can you please point me to any simple pedagogical python examples? Mar 30 at 16:24
• I don't know scipy, so can't help. Mar 30 at 16:49
• You could use sympy to describe the RHS and the Jacobian symbolically, then use the function sympy.lambdify to generate python code that evaluates them and produces numpy arrays, and finally call numpy.linalg.eig. Mar 30 at 19:06