Stiff ODEs

Resources

• LeVeque, Randall J. 2007. “Stiff Ordinary Differential Equations.”

Related

• Ordinary Differential Equations

Main Idea

What does it mean for a system to be stiff?

For an Ordinary Differential Equation $f(y)=y^{\prime }$, linearization $A(y)=\mathrm{\nabla }f(y)$, gives eigenvalues $\lambda$. If the ratio

is large, the system may be stiff (almost a sufficient condition to be stiff, not at all necessary). There's some debate about if this is sufficient as $\lambda$ can have big imaginary components, which means fast oscillation, which requires the numerical method to itself have fine timesteps. However, I think that this is good criteria for many systems. Note that even scalar equations can be stiff. For instance, $y^{\prime }(t)=-y+\sin (500\cdot t)$.

Importantly, systems may be stiff only at some times / states. For instance, a piecewise high frequency forcing term is stiff only when this term is nonzero.

Numerical Methods

For Stiff problems, we want the Absolute Stability Region to extend to large, negative real values, which correspond to a large $max\lambda$ (it's more common to already contain the small $min\lambda$ which is near the origin - zero stability?). Euler's method gives stability region of $\{z\in \mathbb{C}:|z+1|\le 1\}.$ This at most extends to $\mathrm{Re}(z)=-2$, making Euler's method a poor choice for stiff problems. This is an example of a bounded stability region. We want methods that are unbounded (specifically for $z\to -\mathrm{\infty }$). All explicit methods are bounded. Runge-Kutta-Chebychev methods are an exception (still bounded though), and perform better.

Methods that excel at this are stable for the entire left half of the complex plane. These methods are A-stable, and have a stability region that includes at least $\{z\in \mathbb{C}:\mathrm{Re}(z)\le 0\}.$ For instance, backward Euler is A-stable. We can loosen this definition with the definition of A$(\alpha )$-stability. Here, $\alpha$ is the angle above/below the negative real axis where the method is still stable. This is particularly useful if the eigenvalues of our problem have small imaginary components. For instance, A$(0)$-stable is the same as A-stable, and A$(\pi /2)$-stability includes at least the negative real axis.

A (one-step) method is L-stable if it is A-stable and $\underset{z\to \mathrm{\infty }}{lim}|R(z)|=0$. These methods damp transients in effectively one step? We want to use a time step larger than the true decay time of the transient. Backward Euler is L-stable, while Trapezoidal is not. However, Trapezoidal is second-order. Slight changes to the same problem make Backward Euler perform better than Trapezoidal, explained by this L-stability property. This would be an interesting application for Learning Numerics. Note that we can lower the time step to capture the transient more, and then this property does not matter as much.

Runge-Kutta-Chebychev Methods

These are explicit methods that can handle a fair amount of stiffness. Instead of using different stages to increase the order (like the classic 4th order Runge-Kutta, for instance), instead try to increase the stability interval on the negative real axis: maximize $\beta$ , where $[\beta ,0]\subseteq \{z\in \mathbb{C}:R(z)\le 1\}$.

For an $r$-stage RK method, we have the butcher tableau with strictly lower triagonal $A$ and coefficients $b$. These determine $d_i$ in the stability function:

By consistency and first-order (respectively), we need $d_0=1$ and $d_1=1$, which is equivalent to $R(0)=1$ and $R^{\prime }(0)=1$. Next we want to have $|R(\beta )|\le 1$ for maximum $\beta \in \mathbb{R}$. These requirements define a scaled and shifted Chebychev polynomials. Thus, we want our stability function to be the scaled and shifted Chebychev polynomial. This gives $r+2$ values of $z\in \mathbb{R}$ such that $R(z)=1$. Thus, it's easy to add small imaginary perturbations that land us outside of the region of stability. We can add a small margin around the real axis by making $\beta$ smaller. This is done through a damping parameter $w_0>1$. Thus, we've defined the stability function $R(z)$, but still need to come up with the method itself ($A$ and $b$). There are infinitely many solutions. The Runge-Kutta-Chebychev methods specifically use a recurrence relation of the Chebychev polynomials to come up with values. The recurrence relation allows only three function evaluations of $y^{\prime }=f(y)$ to be computed, and the rest of the stages to be constructed from these three. Thus, one can use $r=50$ stages.