11.4. Stiffness

In Section 11.3 we analyzed time step constraints for the semidiscrete heat equation \(\mathbf{u}'=\mathbf{D}_{xx}\mathbf{u}\) in terms of stability regions and the eigenvalues \(\lambda_j\) of the matrix. Since all the eigenvalues are negative and real, the one farthest from the origin, at about \(-4/h^2\), determines the specific time step restriction. For an explicit method, or any method with a finite intersection with the negative real axis, the conclusion is \(\tau=O(h^2)\).

For the Euler and backward Euler solvers, stability is not the only cause of a severely limited time step. Both methods are first-order accurate, with truncation errors that are \(O(\tau)\). Since the spatial discretization we chose is second-order accurate, then we should choose \(\tau=O(h^2)\) for accuracy as well as stability.

This calculus changes for second-order IVP solvers. When both time and space are discretized at second order, the total truncation error is \(O(h^2+\tau^2)\), so it makes sense to use \(\tau=O(h)\) for accuracy reasons alone. Therefore, a stability restriction of \(\tau=O(h^2)\) is much more strict.

Problems for which the time step is dictated by stability rather than accuracy are referred to as stiff. Stiffness is not a binary condition but a spectrum. It can arise in nonlinear problems and in problems having nothing to do with diffusion. Except for the mildest instances of stiffness, an implicit time stepping method is the best choice.


Why should the model equation \(y'=\lambda y\) of absolute stability have wide relevance? Through diagonalization, it is easily generalized to \(\mathbf{u}'=\mathbf{A} \mathbf{u}\) for a constant matrix \(\mathbf{A}\). But that is still a severely limited type of problem.

Consider a general vector nonlinear system


The key to making a connection with absolute stability is to look not at an exact solution but to perturbations of one. Such perturbations always exist in real numerical solutions, such as those due to roundoff error, for example. But if we assume the perturbations are tiny, then we can use linear approximations to describe their evolution. If we conclude from such an approximation that the perturbation may grow without bound, then we must seriously question the value of the numerical solution.

Let’s introduce more precision into the discussion. Suppose that \(\hat{\mathbf{u}}(t)\) is an exact solution that we wish to track, and that a perturbation has pushed us to a nearby solution curve \(\hat{\mathbf{u}}(t) + \mathbf{v}(t)\). Substituting this solution into the governing ODE and appealing to a multidimensional Taylor series, we derive

\[\begin{split}[\hat{\mathbf{u}}(t) + \mathbf{v}(t)]' &= \mathbf{f}\bigl(t,\hat{\mathbf{u}}(t) + \mathbf{v}(t)\bigr), \\ \hat{\mathbf{u}}'(t) + \mathbf{v}'(t) &= \mathbf{f}\left(t, \hat{\mathbf{u}}(t)\right) + \mathbf{J}(t) \mathbf{v}(t) + O\bigl( \|\mathbf{v}(t)\|^2 \bigr).\end{split}\]

We have introduced the Jacobian matrix \(\mathbf{J}\), with entries

(11.4.2)\[J_{ij} = \frac{\partial f_i}{\partial u_j}(t,\hat{\mathbf{u}}(t)).\]

By dropping the higher-order terms, which are negligible at least initially, we derive a linear ODE for the evolution of the perturbation.

Definition 11.4.1 :  Linearization of an ODE

A linearization of system (11.4.1) at an exact solution \(\hat{\mathbf{u}}(t)\) is

(11.4.3)\[\mathbf{v}'(t) = \mathbf{J}(t) \mathbf{v}(t),\]

where \(\mathbf{v}(t)\) is a perturbation to the exact solution, and \(\mathbf{J}\) is the Jacobian matrix (11.4.2).

Example 11.4.2

The Oregonator is a well-known ODE system modeling a chemical oscillator and is given by

(11.4.4)\[\begin{split}\begin{split} u_1' & = s[u_2(1-u_1) + u_1(1-q u_1)], \\ u_2' & = s^{-1}(u_3-u_2-u_1u_2), \\ u_3' & = w(u_1-u_3), \end{split}\end{split}\]

where \(s\), \(q\), and \(w\) are constants. Linearization about an exact (albeit unknown) solution \(\hat{\mathbf{u}}(t)\) leads to the Jacobian

(11.4.5)\[\begin{split}\mathbf{J} (t) = \begin{bmatrix} s(1 -\hat{u}_2 - 2q\hat{u}_1) & s(1-\hat{u}_1) & 0 \\ -\hat{u}_2/s & -(1+\hat{u}_1)/s & 1/s \\ w & 0 & -w \end{bmatrix}.\end{split}\]

Freezing time

While (11.4.3) is linear, the Jacobian matrix in it is time-dependent, which makes analysis difficult. If a perturbation is introduced at a moment \(t=t_\star\), we freeze the Jacobian there and consider

\[\mathbf{v}'=\mathbf{A}\mathbf{v}, \quad \mathbf{A}=\mathbf{J}(t_\star).\]

This equation is of the type we used in Section 11.3 to discuss the absolute stability of IVP solvers. This suggests the following.

Observation 11.4.3 :  Rule of thumb for absolute stability

The eigenvalues of the Jacobian appearing in the linearization about an exact solution, after scaling by the time step \(\tau\), must lie in the stability region of the IVP solver.

We have not stated a theorem here because we made several approximations and assumptions along the way that are not trivial to quantify. Nevertheless, if the rule of thumb is violated, we should expect perturbations to the exact solution to grow significantly with time, eventually rendering the numerical solution useless. Note that roundoff error is constantly introducing perturbations, so the rule of thumb applies along the entire trajectory of the numerical solution.

Demo 11.4.4

In Example 11.4.2 we derived a Jacobian matrix for the Oregonator model. Here is a numerical solution of the ODE.

function ode(u,p,t)
    s,w,q = p
    f = [ 
        s*(u[2]*(1-u[1]) + u[1]*(1-q*u[1])),
        (u[3]-u[2]-u[1]*u[2]) /s,   
    return f

s,w,q = 77.27,.161,8.375e-6
oregon = ODEProblem(ode,[1.,2,3],(0.,500.),[s,w,q])
sol = solve(oregon)
    title="Solution of the Oregonator")
t Solution of the Oregonator

At each value of the numerical solution, we can compute the eigenvalues of the Jacobian. Here we plot all of those eigenvalues in the complex plane.

t,u = sol.t[1:2:end],sol.u[1:2:end]
λ = fill(0.0im,length(t),3)
for (k,u) in enumerate(u)
    J = [
    s*(1-u[2]-2q*u[1]) s*(1-u[1])        0 
         -u[2]/s       -(1+u[1])/s     1/s 
            w               0           -w
    λ[k,:] .= eigvals(J)

    title="Oregonator eigenvalues")
Re(λ) Im(λ) t Oregonator eigenvalues

You can see that there is one eigenvalue that ranges over a wide portion of the negative real axis and dominates stability considerations.

Multiple time scales

The solution to \(y' = \lambda y\), \(y(0)=1\) is \(\exp(\lambda t)\). If \(\lambda\) is real, this solution grows or decays by a factor of \(e\) at \(t=1/|\lambda|\). If \(\lambda = i\omega\) is imaginary, then the solution has sines and cosines of frequency \(\omega\). A complex \(\lambda\) combines these effects.

Observation 11.4.5

We may regard \(|\lambda|^{-1}\), which has units of time, as a characteristic time scale of dynamics due to eigenvalue \(\lambda\).

A Jacobian matrix with eigenvalues at different orders of magnitude therefore implies multiple time scales that the IVP solver needs to cope with. Say \(|\lambda_1|\gg |\lambda_2|\). Any explicit integrator will have a bounded stability region and therefore impose a time step restriction proportional to \(|\lambda_1|^{-1}\). Any good adaptive integrator will obey such a restriction naturally to control the error. But to observe the “slow” part of the solution, the simulation must go on for a time on the order of \(|\lambda_2|^{-1}\), which is much longer.

In Demo 11.4.4, for example, you can see a combination of fast changes and slow evolution.

Demo 11.4.6

The Rodas4P solver is good for stiff problems, and needs few time steps to solve the Oregonator from Demo 11.4.4.

oregon = remake(oregon,tspan=(0.,25.))
sol = solve(oregon,Rodas4P())
println("Number of time steps for Rodas4P: $(length(sol.t)-1)")
Number of time steps for Rodas4P: 127

But if we apply Function 6.5.2 to the problem, the step size will be made small enough to cope with the large negative eigenvalue.

t,u = FNC.rk23(oregon,1e-4)
println("Number of time steps for RK23: $(length(t)-1)")
Number of time steps for RK23: 20773

Starting from the eigenvalues of the Jacobian matrix, we can find an effective \(\zeta(t)\) by multiplying with the local time step size.

λ = fill(1.0im,length(t),3)
for (i,u) in enumerate(u)
    J = [
    s*(1-u[2]-2q*u[1]) s*(1-u[1])        0 
         -u[2]/s       -(1+u[1])/s     1/s 
            w               0           -w
    λ[i,:] .= eigvals(J)

ζ = diff(t).*λ[1:end-1,:]
    title="Oregonator stability")
Re(ζ) Im(ζ) Oregonator stability