Stiff equation

In computational mathematics, a stiff equation is an initial value problem

${\displaystyle {\dot {u}}=f(u)\,,\qquad u(0)=u_{0}\,,\qquad t\in [0,T]\,,}$

where ${\displaystyle f:{\mathbb {R} }^{d}\rightarrow {\mathbb {R} }^{d}}$, requiring dedicated implicit time stepping methods for its efficient numerical integration. The simplest mathematical characterization of a stiff equation is the necessary condition

${\displaystyle {\frac {T}{d}}{\big (}{\mathrm {div} }_{u}\,f{\big )}(u)\ll -1\,.}$

Since ${\displaystyle {\big (}{\mathrm {div} }_{u}\,f{\big )}(u)={\mathrm {trace} }({\mathrm {grad} }_{u}\,f)(u)={\mathrm {trace} }\,f'(u)}$, where ${\displaystyle f'(u)\in {\mathbb {R} }^{d\times d}}$ is the Jacobian matrix of ${\displaystyle f}$ at the point ${\displaystyle u}$, the criterion above is easily evaluated and quantifies stiffness. The criterion is derived, explained and illustrated for nonlinear stiff equations below.

For a linear system with constant coefficients ${\displaystyle {\dot {u}}=Au}$, the divergence is constant, making stiffness a global characteristic whose magnitude is related to the time scale ${\displaystyle T}$.

For a nonlinear system, stiffness usually varies in space and time along the solution trajectory ${\displaystyle u(t)}$, where the criterion quantifies stiffness locally. In practical computations, stiff equations are invariably solved using adaptive methods. ¹²

There is a rich literature on stiff differential equations, but intuitive descriptions and heuristics are far more common than attempts at a rigorous definition of the concept. Hairer and Wanner³ concisely describe the most obvious feature:

Stiff equations are problems for which explicit methods don't work.

This refers to the observation that explicit integration methods are forced to use exceedingly small time steps ${\displaystyle h}$ to maintain numerical stability, preventing such methods from being competitive. Although each step is inexpensive, the total number of steps ${\displaystyle N=T/h}$ becomes prohibitively large, and the integration over ${\displaystyle [0,T]}$ effectively stalls.

By contrast, implicit methods for stiff equations require costly "algebraic" equation solving on each step. The extra work per step is offset by superior stability properties, allowing the use of large time steps. Without crippling stability restrictions, the total computational effort is manageable, and the requested accuracy can be achieved without loss of efficiency. For some stiff equations, the efficiency may be several orders of magnitude higher than that of even the best explicit methods.

The phenomenon of stiffness is analogous to well-known performance issues in other areas of numerical analysis. For example, in optimization, gradient methods have similar limitations. Using the steepest descent method to minimize a convex functional ${\displaystyle F(x)}$ leads to the iteration

${\displaystyle x^{k+1}=x^{k}-h\,{\mathrm {grad} }_{x}F(x^{k})\,,}$

where ${\displaystyle h}$ is the line search step size. This is merely the explicit Euler method for the differential equation ${\displaystyle {\dot {x}}=-{\mathrm {grad} }_{x}F(x^{k})}$. Problems with steep gradients are known for their slow convergence, due to restrictions on the choice of ${\displaystyle h}$. The typical remedy is to replace the gradient method by some Newton-type method, which essentially corresponds to using an "implicit" method for differential equations.

Another example is in solving nonlinear equations of the form ${\displaystyle x=g(x)}$. Simple fixed point iteration, ${\displaystyle x^{k+1}=g(x^{k})}$, does not converge for problems where the Lipschitz constant ${\displaystyle L[g]}$ is large, but only for contractions. Yet the equation may have a unique solution under much weaker but common conditions, e.g. when the map ${\displaystyle g-I}$ is monotone and the logarithmic Lipschitz constant of ${\displaystyle g}$ satisfies ${\displaystyle M[g]<1}$. Since fixed point iterations diverge, it is necessary to use a Newton-type iteration.

A third example is iterative solvers for the (discretized) Poisson equation ${\displaystyle -\Delta u=f}$. Here, the Jacobi method (which attempts to avoid matrix algebra) is equivalent to using the explicit Euler method (in pseudo time) for the corresponding method-of-lines discretization of the diffusion equation ${\displaystyle u_{t}=\Delta u+f}$, seeking its stationary solution. The diffusion equation is a prototypical example of a stiff equation: the time step of the explicit recursion is severely restricted by a CFL condition and required to be proportional to the square of the mesh width in space. This makes the convergence of the Jacobi iteration painfully slow. As before the remedy is to use some form of algebraic equation solving, e.g. by using a preconditioner or a multigrid method.

Thus there are many different mathematical problems, where computational efficiency necessitates the use of dedicated numerical methods, typically of a considerable complexity. Stiff equations is the particular instance for initial value problems in ordinary differential equations. Nevertheless, with a well designed adaptive stiff solver, the efficient numerical integration of most stiff equations is now a routine task.

The first mention of stiff equations is found in Curtiss and Hirschfelder in 1952 ⁴. The authors discuss a scalar model equation equivalent to

${\displaystyle {\dot {x}}={\frac {1}{a(t,x)}}{\big (}x-g(t){\big )}}$

and characterize stiffness as follows:

If ${\displaystyle \Delta t}$ is the desired resolution of ${\displaystyle t}$ or the interval which will be used in the numerical integration, the equation is "stiff" if

${\displaystyle \left|{\frac {a(t,x)}{\Delta t}}\right|\ll 1}$

and ${\displaystyle g}$ is well behaved.

This characterization relates a time scale ${\displaystyle \Delta t}$ to the decay rate of transients, governed by the coefficient ${\displaystyle a(t,x)}$, which is assumed to be negative. As ${\displaystyle a(t,x)}$ varies in time and space, it allows stiffness to vary accordingly. Prothero and Robinson introduced a similar but more instructive model problem,⁵ essentially taking ${\displaystyle a(t,x)=1/\lambda <0}$ constant and studying the equation

${\displaystyle {\dot {x}}=\lambda {\big (}x-g(t){\big )}+{\dot {g}}(t)\,,\qquad x(0)=x_{0}\neq g(0)\,,\qquad t\in [0,T].}$

The solution consists of a decaying transient ${\displaystyle {\big (}x_{0}-g(0){\big )}\,{\mathrm {e} }^{t\lambda }}$, together with the particular solution ${\displaystyle g(t)}$, which is supposed to be bounded and smooth. The mathematical solution ${\displaystyle x(t)}$ is then obviously bounded. Contrasting the explicit and implicit Euler methods, one can investigate under what circumstances the numerical methods replicate the properties of the mathematical solution. The Curtiss-Hirschfelder criterion will be considered for ${\displaystyle \Delta t=T}$ and ${\displaystyle \Delta t=h}$, corresponding to the conditions ${\displaystyle T\lambda \ll -1}$ and ${\displaystyle h\lambda \ll -1}$, respectively.

The explicit Euler method for a differential equation ${\displaystyle {\dot {u}}=f(u)}$ is defined by

${\displaystyle u_{n+1}=u_{n}+hf(u_{n})\,,}$

while the implicit Euler method is defined by the recursion

${\displaystyle u_{n+1}=u_{n}+hf(u_{n+1})\,.}$

Here ${\displaystyle h}$ is the (constant) time step, and ${\displaystyle u_{n}}$ approximates the exact solution ${\displaystyle u(t_{n})}$ at time ${\displaystyle t_{n}=nh}$. In the explicit method ${\displaystyle u_{n+1}}$ is computed directly by an evaluation of the vector field ${\displaystyle f}$. In the implicit method, however, one must solve the "algebraic" equation ${\displaystyle u_{n+1}-hf(u_{n+1})=u_{n}}$ for ${\displaystyle u_{n+1}}$. This improves the stability of the recursion.

For the Prothero-Robinson problem the explicit Euler method produces the numerical solution (on integrated form)

${\displaystyle x_{n}=\left(1+h\lambda \right)^{n}x_{0}+\sum _{k=0}^{n-1}\left(1+h\lambda \right)^{n-k-1}h{\big (}{\dot {g}}(t_{k})-\lambda g(t_{k}){\big )}\,.}$

The impediment is clearly visible: positive powers of the factor ${\displaystyle 1+h\lambda }$ produce exponential growth (numerical instability) unless we impose the stability condition ${\displaystyle |1+h\lambda |\leq 1}$. Since ${\displaystyle \lambda }$ is real and negative, this implies ${\displaystyle -2\leq h\lambda \leq 0}$. Thus the explicit Euler method cannot use step sizes adapted to the regularity of ${\displaystyle g(t)}$. Instead the total number of steps to complete the integration, ${\displaystyle N=T/h\sim T|\lambda |}$, becomes extremely large. Efficiency is inversely proportional to stiffness as quantified by the product ${\displaystyle T|\lambda |\gg 1}$.

As a consequence the explicit Euler method may take "forever" to solve a scalar problem with a perfectly smooth mathematical solution ${\displaystyle \sim g(t)}$ once the (fast) transient has decayed. This is also the case even if the initial value is ${\displaystyle x(0)=g(0)}$. Thus, the step size restriction is not a matter of whether the transient is present. It is due to the fact that any step size ${\displaystyle h}$ such that ${\displaystyle |1+h\lambda |>1}$ inevitably excites numerical instability.

If instead the implicit Euler method is used, the numerical solution becomes

${\displaystyle x_{n}=\left(1-h\lambda \right)^{-n}x_{0}+\sum _{k=0}^{n-1}\left(1-h\lambda \right)^{k-n}h{\big (}{\dot {g}}(t_{k+1})-\lambda g(t_{k+1}){\big )}\,.}$

Here, the stability requirement is ${\displaystyle |1-h\lambda |^{-1}\leq 1}$, which is satisfied for all ${\displaystyle h>0}$ when ${\displaystyle {\mathrm {Re} }\,\lambda <0}$, no matter how large ${\displaystyle h|\lambda |}$ is. A minor issue is that in practice one needs to resolve the transient (this may initially call for a short step size), but once it has decayed, one can use "large" step sizes, in the sense that ${\displaystyle h}$ is adapted to the regularity of ${\displaystyle g(t)}$, irrespective of the magnitude of ${\displaystyle |h\lambda |}$. Of the two methods above, only the implicit Euler method can handle stiff equations efficiently.

Example Consider the Prothero-Robinson test problem

${\displaystyle {\dot {x}}\,=\,\lambda (x-\sin \omega t)+\omega \cos \omega t\,;\qquad x(0)=x_{0}\,,}$

with exact solution ${\displaystyle x(t)=x_{0}\,{\mathrm {e} }^{\lambda t}+\sin \omega t}$. Taking ${\displaystyle x_{0}=0}$, the homogeneous solution is not present in the exact solution ${\displaystyle x(t)=\sin \omega t}$, but the exponential term will show up in the local solutions passing through the points ${\displaystyle x_{n}}$ generated by the numerical methods. Taking ${\displaystyle \omega =\pi }$ the problem is solved on ${\displaystyle [0,1]}$, using ${\displaystyle N=5}$ and ${\displaystyle 10}$ steps, respectively (${\displaystyle h=0.2}$ and ${\displaystyle h=0.1}$). The experiment demonstrates the impact of varying ${\displaystyle h\lambda }$, which affects the stability/damping of the numerical methods. We take ${\displaystyle \lambda =-20}$, using the same computational setup for the explicit and implicit Euler methods.

The results are shown in Figure 1. With ${\displaystyle T\lambda =-20}$, there is a moderately strong exponential damping, as is evident from the fast contracting local solutions. The setup is deliberately understated for visualization purposes; in real stiff equations much larger values would be typical. Nevertheless, ${\displaystyle \lambda }$ is sufficiently large and negative to demonstrate the step size limitations of the explicit Euler method. For ${\displaystyle N=5}$, its computed solution becomes oscillatory and diverges from the exact solution, even though the initial value was taken on the exact solution. Repeating the same exercise for ${\displaystyle T\lambda =-2000}$, say, shows how the difference between the two methods grows with increasing stiffness.

In more realistic problems, the instability is often more dramatic and referred to as the solution "blowing up". Here, the maximum stable step size (${\displaystyle h=0.1}$) has only been exceeded by a small factor. Meanwhile, the implicit Euler method proceeds without loss of stability, and is able to produce accurate results also for the larger step size.

In real computations, stiff equations are solved using adaptive methods. These select the step size automatically to keep the computational process in check, usually controlling both accuracy and stability. Thus, if an adaptive explicit method is used, the solution may not go unstable, as the adaptivity takes measures not to exceed the maximum stable step size. Instead, the method proceeds at the (prohibitive) cost of low efficiency.

The essence of the analysis is that a problem parameter, ${\displaystyle \lambda }$, which is characteristic of the vector field and related to the transient behavior of solutions, is linked to a time scale, either the step size ${\displaystyle h}$ or the range of integration ${\displaystyle T}$ directly. The products ${\displaystyle h\lambda }$ or ${\displaystyle T\lambda }$ are both dimensionless, scaling invariant, and negative. A general stiffness concept must reflect these properties, and also be able to address stiffness in nonlinear differential equations.

The scalar Prothero-Robinson problem can be extended to linear systems of equations,

${\displaystyle {\dot {x}}=A{\big (}x-g(t){\big )}+{\dot {g}}(t)\,,\quad A\in {\mathbb {R} }^{d\times d}}$

as well as nonlinear systems,

${\displaystyle {\dot {x}}=f{\big (}x-g(t){\big )}+{\dot {g}}(t)\,,\quad f:{\mathbb {R} }^{d}\rightarrow {\mathbb {R} }^{d}}$

provided that ${\displaystyle f(0)=0}$. The general stiffness concept needs to specify what is meant by a vector field being "large and negative". The analysis above also shows that the smooth function ${\displaystyle g}$ at most plays a minor role. Thus, by substituting ${\displaystyle u}$ for ${\displaystyle x-g(t)}$, we have three types of problems,

${\displaystyle {\begin{aligned}{\dot {u}}\,&=\,\lambda u\\{\dot {u}}\,&=\,Au\\{\dot {u}}\,&=\,f(u)\,.\end{aligned}}}$

Here stiffness needs to be characterized for matrices ${\displaystyle A}$ as well as for nonlinear maps ${\displaystyle f}$. The first equation above is the well known Dahlquist linear test equation, primarily used to determine the stability region of a time stepping method. The stability region ${\displaystyle S}$ of a method is the set of all ${\displaystyle h\lambda \in \mathbb {C} }$ such that the numerical method produces bounded solutions. The previous discussion of stiffness is effectively based on the Dahlquist test equation.

In order to consider the second, linear system of equations, it is common to make use of the eigenvalues ${\displaystyle \lambda _{k}}$ of the matrix ${\displaystyle A}$. After all, for a numerical method to produce bounded (stable) solutions to ${\displaystyle {\dot {u}}=Au}$, it is necessary and sufficient that all eigenvalues satisfy ${\displaystyle h\lambda _{k}\in S}$. The previous discussion carries over by considering one eigenvalue at a time. The linear problem is then stiff if all eigenvalues have negative real parts, and at least one of them satisfies ${\displaystyle T\lambda _{k}\ll -1}$; the latter eigenvalue will then cause the typical step size restrictions experienced by an explicit method, due to the fact that all explicit methods have bounded stability regions, and require ${\displaystyle h\lambda _{k}\in S}$. The step size has to be selected small enough to fulfill this requirement, see Figure 2.

Yet, it is common to find a different characterization of stiffness, in terms of the "stiffness ratio", defined by

${\displaystyle {\frac {\max _{k}|{\mathrm {Re} }\,\lambda _{k}|}{\min _{k}|{\mathrm {Re} }\,\lambda _{k}|}}\,,}$

assuming that ${\displaystyle \lambda _{k}\in {\mathbb {C} }^{-}}$.⁶ Although stiff systems often do have a large stiffness ratio, Byrne and Hindmarsh emphasize that this is neither necessary nor sufficient for ${\displaystyle {\dot {u}}=Au}$ to be a stiff equation.⁷

First, there is no comparison with the time scale ${\displaystyle T}$, but only an intrinsic comparison of "time constants" ${\displaystyle 1/\lambda _{k}}$, where a large stiffness ratio merely implies that there are "faster" and "slower" transient components, often referred to as widely different time scales. Second, as noted above, there are scalar stiff equations; these have stiffness ratio ${\displaystyle 1}$. And third, the stiffness ratio breaks down for problems with a singular matrix ${\displaystyle A}$, which is common in some applications. (An example is chemical reaction kinetics, where the singularity corresponds to conservation of mass.)

Thus the stiffness ratio is misleading: it fails to characterize the most important aspects of stiff equations.

In a similar vein, it has been suggested that the nonlinear problem is stiff if the Jacobian matrix ${\displaystyle f'(u)}$ has a large stiffness ratio.⁸ Naturally, this too fails, as noted by Artemiev and Averina:⁹

For example, a nonlinear autonomous ODE system with ${\displaystyle f(0)=0}$ can be called stiff if real parts [sic] of the eigenvalues of the Jacobi matrix [${\displaystyle f'(0)}$] satisfy the above conditions. However, the famous van der Pol equation, which is frequently used as an example of stiff ODEs, does not satisfy this definition ... it is impossible to determine stiffness only by means of eigenvalues of the Jacobi matrix for nonlinear systems

Although in part valid, this argument too is off the mark, as the zero solution of the van der Pol equation is an unstable equilibrium. Thus stiffness does not occur near the origin, but only along the limit cycle as will be seen below, where the van der Pol equation is analyzed in full detail. A precise stiffness characterization must therefore be local, capturing the behavior of (small) perturbations near the actual solution.

Replacing the stiffness ratio, the modern approach is to construct a functional ${\displaystyle s_{2}[\cdot ]}$, called the stiffness indicator, such that the criteria ${\displaystyle s_{2}[A]\!\cdot \!T\ll -1}$ and ${\displaystyle s_{2}[f]\!\cdot \!T\ll -1}$, respectively, quantify stiffness by relating a distinct vector field property to the time scale ${\displaystyle T}$. This characterization is based on additional observations. Thus Dekker and Verwer¹⁰ write (original emphasis):

The essence of stiffness is that the solution to be computed is slowly varying but that perturbations exist which are rapidly damped.

Shampine¹¹ expresses a similar view, but is conceptually more specific, relating stiffness to dissipation, and what in practice is irreversibility:

A way we prefer for describing the latter condition is that [the solution] is very unstable in the opposite direction [of time].

The reverse instability can easily be characterized using logarithmic norms. Thus, for ${\displaystyle {\dot {u}}=Au}$, the norm of the solution can be bounded using differential inequalities. Let ${\displaystyle \|u\|_{2}^{2}=u^{*}u}$ be the Euclidean norm, possibly with ${\displaystyle u\in {\mathbb {C} }^{d}}$. Then

${\displaystyle m_{2}[A]\cdot \|u\|_{2}\,\leq \,{\mathrm {D} }_{t}\,\|u\|_{2}\,\leq \,M_{2}[A]\cdot \|u\|_{2}\,,}$

where ${\displaystyle {\mathrm {D} }_{t}}$ denotes the time derivative and ${\displaystyle m_{2}[A]}$ and ${\displaystyle M_{2}[A]}$ are the lower and upper logarithmic norms of ${\displaystyle A}$, respectively. These are the extrema of a symmetric quadratic form,

${\displaystyle m_{2}[A]\,\leq \,{\frac {u^{*}{\mathrm {He} }(A)u}{u^{*}u}}\,\leq \,M_{2}[A]\,,}$

where ${\displaystyle {\mathrm {He} }(A)=(A+A^{*})/2}$ denotes the Hermitian part of ${\displaystyle A}$. Its stationary points are found from the symmetric eigenvalue problem

${\displaystyle {\mathrm {He} }(A)v=\mu v\,,}$

whose ${\displaystyle d}$ real eigenvalues ${\displaystyle \mu _{1}\geq \mu _{2}\geq \dots \geq \mu _{d}}$ are referred to as the logarithmic values of ${\displaystyle A}$. In relation to the eigenvalues ${\displaystyle \lambda _{k}}$ of ${\displaystyle A}$, the logarithmic values satisfy

${\displaystyle m_{2}[A]\,\equiv \,\mu _{d}\,\leq \,\min _{k}{\mathrm {Re} }\,\lambda _{k}\,;\qquad \max _{k}{\mathrm {Re} }\,\lambda _{k}\,\leq \,\mu _{1}\,\equiv \,M_{2}[A]\,,}$

together with the trace identity

${\displaystyle {\mathrm {Re} }\,{\mathrm {trace} }\,A\,=\,\sum _{k}{\mathrm {Re} }\,\lambda _{k}\,=\,\sum _{k}\mu _{k}\,=\,{\mathrm {trace} }\,{\mathrm {He} }(A)\,.}$

Since ${\displaystyle u(t)={\mathrm {e} }^{tA}u_{0}}$ it follows from the differential inequality above that

${\displaystyle {\mathrm {e} }^{tm_{2}[A]}\,\leq \,\|{\mathrm {e} }^{tA}\|_{2}\,\leq \,{\mathrm {e} }^{tM_{2}[A]}\,\qquad t\geq 0\,,}$

bounding the maximum growth and decay rates in the system. Thus, over a time interval of length ${\displaystyle T>0}$, the maximum possible forward time growth is ${\displaystyle {\mathrm {e} }^{TM_{2}[A]}}$, while the maximum reverse time growth is ${\displaystyle {\mathrm {e} }^{-Tm_{2}[A]}}$.

Söderlind et al.¹² now quantify Shampine's observation. Thus, in a stiff system,

${\displaystyle {\mathrm {e} }^{-Tm_{2}[A]}\,\gg \,{\mathrm {e} }^{TM_{2}[A]}\,,}$

from which it follows that ${\displaystyle {\mathrm {e} }^{T{\big (}m_{2}[A]+M_{2}[A]{\big )}}\ll 1}$, i.e., ${\displaystyle T{\big (}m_{2}[A]+M_{2}[A]{\big )}\ll -1}$. The stiffness indicator is defined¹³

${\displaystyle s_{2}[A]\,=\,{\frac {m_{2}[A]+M_{2}[A]}{2}}\,,}$

and a stiff equation is characterized by the condition

${\displaystyle s_{2}[A]\cdot T\ll -1\,.}$

The reason for including ${\displaystyle M_{2}[A]}$ in the definition of the stiffness indicator is that it may often happen in nonlinear systems that ${\displaystyle M_{2}[A]}$ is positive, offsetting stiffness. This is then accounted for by defining the stiffness indicator so as to have odd parity, i.e., ${\displaystyle s_{2}[-A]=-s_{2}[A]}$.

Shampine's criterion implies that ${\displaystyle m_{2}[A]}$ is large and negative, i.e., the lower bound on ${\displaystyle \|{\mathrm {e} }^{tA}\|_{2}}$ is extremely small. In other words, the flow is nearly a semigroup, and ${\displaystyle {\mathrm {e} }^{tA}}$ is "close to singular". Unlike in a Lipschitz system, where one can solve (theoretically) the differential equation also in reverse time, this is effectively impossible in a stiff system. The reverse time stiff equation corresponds to an extremely ill conditioned "inverse problem".

Using logarithmic norms for nonlinear maps¹⁴ the same arguments apply to a nonlinear vector field ${\displaystyle f}$, but since ${\displaystyle s_{2}[f]}$ is defined globally, the stiffness indicator is instead defined locally as ${\displaystyle s_{2}[f'(u)]}$ along the trajectory.

In addition to the stiffness indicator, Söderlind et al.¹⁵ introduce a local reference time scale ${\displaystyle \Delta t(u)}$ defined by

${\displaystyle \Delta t(u)\,=\,1/\max {\big (}1/T,-s_{2}[f'(u)]{\big )}\,.}$

Stiffness is then characterized by the method independent condition ${\displaystyle T/\Delta t(u)\gg 1}$, referred to as the (local) stiffness factor. A large stiffness factor is a necessary condition for stiffness. But method dependent criteria can also be assessed. In numerical computations, ${\displaystyle \Delta t(u)}$ can be compared to the actual step size ${\displaystyle h}$, whose local step size stiffness factor ${\displaystyle h/\Delta t(u)}$ depends on the choice of integration method and the accuracy requirement. Only implicit methods can use step sizes with ${\displaystyle h/\Delta t(u)\gg 1}$, but there may also be parts of the integration where the step size stiffness factor is only moderate.

An important aspect of implicit time stepping methods is the need to solve equations of the form

${\displaystyle u=\gamma hf(u)+\psi }$

for ${\displaystyle u}$ on every step. Here ${\displaystyle \gamma >0}$ is a method characteristic constant of moderate size, and ${\displaystyle \psi }$ is a known vector. The question arises whether this equation has a (unique) solution when ${\displaystyle hf}$ is "large", in the sense that ${\displaystyle hs_{2}[f']\ll -1}$, which represents stiff computation. Rewriting the equation as

${\displaystyle {\big (}\gamma hf-I{\big )}(u)=-\psi \,,}$

the uniform monotonicity theorem applies: the equation has a unique solution if the map ${\displaystyle \gamma hf-I}$ is (negative) monotone, i.e., if ${\displaystyle M_{2}[\gamma hf-I]<0}$. Since

${\displaystyle M_{2}[\gamma hf-I]<0\quad \Leftrightarrow \quad M_{2}[\gamma hf]<1\,,}$

this condition is usually satisfied by a fair margin in stiff computation. If ${\displaystyle M_{2}[f]<0}$, it is satisfied for all ${\displaystyle h>0}$, but if ${\displaystyle M_{2}[f]>0}$, it may require a minor step size restriction. In the Shampine criterion, as well as in the stiffness indicator, it holds that

${\displaystyle {\frac {m_{2}[\gamma hf]+M_{2}[\gamma hf]}{2}}\,=\,\gamma h\!\cdot \!s_{2}[f]\ll -1\,,}$

where ${\displaystyle m_{2}[\gamma hf]\ll -1}$, and where ${\displaystyle M_{2}[\gamma hf]}$ is moderate if it is positive. Therefore, although the equation typically requires Newton iteration in stiff computation (${\displaystyle \gamma hf}$ is not a contraction), existence and uniqueness are usually guaranteed also when the step size stiffness factor is large.

Because it is relatively expensive to compute the stiffness indicator by solving the symmetric eigenvalue problem ${\displaystyle \,{\mathrm {He} }{\big (}f'(u){\big )}v=\mu v\,}$ for ${\displaystyle \mu _{1}}$ and ${\displaystyle \mu _{d}}$, an inexpensive alternative is of interest. Noting that the stiffness indicator is the arithmetic average of the largest and the smallest logarithmic value, an option is to replace this average by the average of all logarithmic values. With the aforementioned trace identity (see above), we define the scaled (arithmetic average) real trace functional¹⁶

${\displaystyle \tau [A]\,=\,{\frac {1}{d}}\,{\mathrm {Re} }\,{\mathrm {trace} }(A)\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}{\mathrm {Re} }\,a_{kk}\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}{\mathrm {Re} }\,\lambda _{k}\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}\mu _{k}\,,}$

noting that the computation of ${\displaystyle \tau [A]}$ is both inexpensive and trivial: no eigenvalue computations are required. The scaled trace (as well as the divergence) is a linear functional, hence ${\displaystyle \tau [-A]=-\tau [A]}$, sharing the odd parity of ${\displaystyle s_{2}[A]}$. Moreover, for a nonlinear map,

${\displaystyle \tau [f'(u)]\,=\,{\frac {1}{d}}\,{\mathrm {trace} }\,f'(u)\,=\,{\frac {1}{d}}\,{\mathrm {trace} }({\mathrm {grad} }_{u}\,f)(u)\,=\,{\frac {1}{d}}\,({\mathrm {div} }_{u}\,f)(u)\,.}$

The functional ${\displaystyle \,\tau [f'(u)]\,}$ serves as an alternative stiffness indicator. Stiffness can therefore also be characterized by the divergence condition

${\displaystyle {\frac {T}{d}}\,{\big (}{\mathrm {div} }_{u}\,f{\big )}(u)\,\ll \,-1\,,}$

which is readily evaluated for most problems.

For a linear system ${\displaystyle {\dot {u}}=Au}$, it holds that

${\displaystyle {\mathrm {D} }_{t}\,{\mathrm {det} }\,{\mathrm {e} }^{tA}\,=\,{\mathrm {trace} }\,A\cdot {\mathrm {det} }\,{\mathrm {e} }^{tA}\,.}$

Defining a scaled (geometric average) absolute determinant ${\displaystyle \,\delta [A]=|{\mathrm {det} }\,A|^{1/d}\,}$ we have ${\displaystyle \,\delta [I]=1\,}$ and ${\displaystyle \,\delta [\alpha A]=|\alpha |\!\cdot \!\delta [A]\,}$, avoiding the standard determinant's direct dependence on the dimension ${\displaystyle d}$. (This is of special importance when the differential equation is derived from a method-of-lines discretization of a partial differential equation.) Here ${\displaystyle \,\delta [A]\,}$ is the linear "size" of the phase space volume spanned by the column vectors of ${\displaystyle A}$, such as in ${\displaystyle \,\delta [2I]=2\,}$, independent of ${\displaystyle d}$. The differential equation for the determinant of the flow can now be rewritten

${\displaystyle {\mathrm {D} }_{t}\,\delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,=\,\tau [A]\cdot \delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,.}$

Hence

${\displaystyle \delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,=\,{\mathrm {e} }^{t\,\tau [A]}\,,}$

expressing that an initial (linear) phase volume changes by a factor ${\displaystyle \,\delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,}$ over a time ${\displaystyle t}$. If ${\displaystyle \,\tau [A]<0\,}$, phase volume is compressed by the flow ${\displaystyle \,{\mathrm {e} }^{tA}\,}$ for ${\displaystyle t\geq 0}$.

It follows that ${\displaystyle \,\tau [A]\leq 0\,}$ is a necessary condition for stability. (Compare the sufficient criterion ${\displaystyle \,M_{2}[A]\leq 0\,}$ for the stability of the zero solution.) Shampine's characterization of stiffness, in terms of strongly unstable solutions in reverse time, can be written as ${\displaystyle \,T\!\cdot \tau [-A]\gg 1\,}$, which, due to the odd parity of the trace, is equivalent to the stiffness condition ${\displaystyle T\!\cdot \tau [A]\ll -1}$. These observations, concepts and constructions provide the key to characterizing and quantifying stiffness, with the divergence criterion ${\displaystyle \,T\!\cdot \tau [f']\ll -1\,}$ being the simplest.

A large number of nontrivial test problems have been collected in the Bari Test Set for Initial Value Problem Solvers.¹⁷ This contains detailed examples of problems, solvers, and results achieved under specified conditions by well-known codes.

Here three stiff nonlinear problems are selected to illustrate how the (simplified) stiffness indicator is used to identify and characterize stiffness correctly. This includes the ability to distinguish the complex and rapidly changing behavior at (locally unstable) turning points, such as those observed in the van der Pol and Oregonator equations. Further, the stiffness indicator also distinguishes between subintervals where stiffness is prominent and where the problem is nonstiff. Finally in some cases it is even possible to identify which of the dependent variables contributes most to stiffness.

Operational criteria, such as how the adaptive time step depends on the accuracy requirements and varies along the solution, are also illustrated while quantifying the step size stiffness factor accordingly. Taken together, a detailed insight can be derived from the stiffness indicator, some of it even before the actual computation begins.

Problem 1. (Flame propagation) Consider the scalar flame propagation model¹⁸ for ${\displaystyle t\in [0,1]}$,

${\displaystyle \varepsilon \,{\dot {u}}\,=\,2\,(u^{2}-u^{3})\,;\qquad u(0)=\varepsilon \ll 1\,.}$

This model is commonly used as a test problem to challenge the time step adaptivity of software for stiff equations.

For ${\displaystyle u\in (0,1)}$, the right hand side is positive and the solution ${\displaystyle u(t)}$ is monotonically increasing. Starting at ${\displaystyle u(0)=\varepsilon }$, it moves away from the unstable equilibrium at ${\displaystyle u=0}$. The divergence of the vector field (stiffness indicator) is simply ${\displaystyle f'(u)=2\,(2u-3u^{2})/\varepsilon }$. Since

${\displaystyle f'(\epsilon )\approx 4\,;\qquad f'(1)=-\,{\frac {2}{\varepsilon }}\ll -1\,,}$

the solution starts out nonstiff, when ${\displaystyle u\ll 1}$. At ${\displaystyle t\approx 0.5}$ when ${\displaystyle u\sim 1/2}$, instability has become catastrophic, with an extremely rapid transition towards the second, stable equilibrium ${\displaystyle u=1}$.

When ${\displaystyle u>2/3}$, stability is regained, and the solution enters the stiff regime, see Figure 3, where we have used ${\displaystyle \varepsilon =10^{-3}}$. The equation is stiffer the smaller ${\displaystyle \,\varepsilon \,}$ is chosen, and the closer to ${\displaystyle \,u=1\,}$ the solution moves.

The demonstration compares an explicit third order Runge-Kutta method with a second order error estimator, with a third order A-stable implicit Runge-Kutta method, also with a second order error estimator. The implicit method has no advantage in the nonstiff regime, but in the stiff regime it can use much larger steps. The difference between the two methods increases for smaller ${\displaystyle \,\varepsilon }$.

Problem 2. (van der Pol equation) The van der Pol equation is usually presented in a couple of alternate forms, with different scalings. It is a nonlinear problem whose solutions approach a limit cycle. Here we scale time to write the system as

${\displaystyle {\begin{aligned}{\dot {x}}&=2\kappa \cdot y\\{\dot {y}}&=2\kappa ^{2}\cdot (1-x^{2})\,y-2\kappa \cdot x\,,\end{aligned}}}$

with initial conditions ${\displaystyle x(0)=2}$, ${\displaystyle y(0)=0}$ chosen on the limit cycle. The time scaling has been chosen so as to make the period of the limit cycle ${\displaystyle \mathrm {O} (1)}$ nearly independent of the parameter ${\displaystyle \kappa }$. The problem can then be solved on ${\displaystyle [0,1]}$, accounting for a full period. The problem originates in the study of nonlinear oscillations in electric circuits.

For large values of ${\displaystyle \kappa }$ (see Figure 4, where ${\displaystyle \kappa =200}$), the limit cycle consists of two stable branches, connected by nearly discontinuous transitions showing a catastrophic loss of stability. Due to the trace identity, the stiffness indicator ${\displaystyle s_{2}[f'(u)]}$ coincides with the scaled trace ${\displaystyle \tau [f'(u)]}$, given by

${\displaystyle s_{2}[f'(u)]\,=\,\tau [f'(u)]\,=\,\kappa ^{2}\,(1-x^{2})\,,}$

where ${\displaystyle u}$ is the vector ${\displaystyle (x,y)^{\mathrm {T} }}$. The stiffness indicator is readily identified in the second equation of the van der Pol system, hiding in plain sight.

The scaled trace has a constant term, ${\displaystyle \kappa ^{2}}$, which represents the contribution from the linear part of the system. As this term is positive it is destabilizing. The equilibrium ${\displaystyle x=y=0}$ is unstable, since ${\displaystyle \tau [f'(0)]=\kappa ^{2}>0}$. The scaled trace also has a negative term ${\displaystyle -\kappa ^{2}x^{2}}$ representing the nonlinear contribution to stiffness. As the scaled trace is independent of ${\displaystyle y}$, the latter variable has no particular influence on stiffness.

The solution enters the stiff regime when ${\displaystyle \tau [f'(u)]\ll -1}$. Since ${\displaystyle \tau [f'(u)]=\kappa ^{2}\,(1-x^{2})}$ it follows that stiffness only occurs for ${\displaystyle \,|x|>1}$, i.e., along the two stable branches, where stiffness is proportional to ${\displaystyle \,\kappa ^{2}}$, and can be arbitrarily large. In particular, it follows that if an explicit method were to be used to solve the problem along the stable branch, the computational effort is ${\displaystyle \mathrm {O} (\kappa ^{2})}$. By contrast, the effort is ${\displaystyle \mathrm {O} (1)}$ for a stiff solver.¹⁹

Note that all this information is available a priori. It is confirmed by the numerical solution of the problem.

Problem 3. (Oregonator equation) Chemical reaction kinetics is a rich source of stiff differential equations. The Oregonator equation models an autocatalytic reaction involving three reactions. For suitable initial conditions, it has a limit cycle of period ${\displaystyle T\approx 305}$. The equations are

${\displaystyle {\begin{aligned}{\dot {x}}\,&=\,s\cdot (x-xy+y-q\,x^{2})\\{\dot {y}}\,&=\,(z-y-xy)/s\\{\dot {z}}\,&=\,w\cdot (x-z)\,.\end{aligned}}}$

The test problem uses the parameters ${\displaystyle s=77.27}$, ${\displaystyle q=8.375\cdot 10^{-6}}$, ${\displaystyle w=0.161}$, together with initial conditions ${\displaystyle x(0)=3.2,\,y(0)=1.45,\,z(0)=2.5}$.²⁰

Because the dimension of the system is ${\displaystyle d=3}$, the stiffness indicator ${\displaystyle s_{2}[f']}$ and the scaled trace (divergence) ${\displaystyle \tau [f'(u)}$ do not coincide, although the differences are small. More importantly, because the equation is a more intricate example of stiffness, we need to gain preliminary insight by keeping the parameters as long as possible in the analysis. We therefore choose the scaled divergence, which offers this advantage. Thus we find

${\displaystyle T\cdot \tau [f'(u)]\,=\,\alpha +\beta \,x+\gamma \,y,}$

where

${\displaystyle {\begin{aligned}\alpha \,&=\,T\cdot (s^{2}-sw-1)/(3s)\,\approx \,7.838\cdot 10^{3}\\\beta \,&=\,-T\cdot (2qs^{2}+1)/(3s)\,\approx \,-1.447\\\gamma \,&=\,-T\cdot s/3\,\approx \,-7.856\cdot 10^{3}.\end{aligned}}}$

Here the constant term ${\displaystyle \alpha }$ is the linear contribution to the divergence. It is positive, therefore destabilizing. The two other terms represent nonlinear contributions, as they multiply ${\displaystyle x}$ and ${\displaystyle y}$. Their coefficients ${\displaystyle \beta }$ and ${\displaystyle \gamma }$ are negative, and will contribute to stiffness. But as the divergence is independent of ${\displaystyle z}$, the latter variable does not contribute to stiffness other than to the constant term ${\displaystyle \alpha }$. The largest coefficient is ${\displaystyle \gamma }$, indicating that stiffness will be particularly strong when ${\displaystyle y}$ is large and positive. Whether ${\displaystyle x}$ will contribute depends on what magnitude it will reach. In Figure 5, the actual data are collected during one full period ${\displaystyle [0,305]}$ plotted on the normalized interval ${\displaystyle [0,1]}$.

During a brief interval when ${\displaystyle x}$ is large, the scaled divergence is of order ${\displaystyle T\cdot \tau [f'(u)]\approx -10^{5}}$ (just barely visible as a notch in the upper left corner of the high resolution scaled trace graph), making the problem stiff there. But the worse is yet to come; when ${\displaystyle y}$ becomes large stiffness becomes severe, with scaled divergence at ${\displaystyle T\cdot \tau [f'(u)]\approx -10^{7}}$.

When an adaptive third order A-stable Runge-Kutta method is used, the step sizes exceed the reference time scale by three to four orders of magnitude. This corresponds to the efficiency gain of the stiff solver over the potential use of an explicit method.

As noted above, stiff equations require implicit time stepping methods, either Runge-Kutta methods or linear multistep methods. There are also alternatives based on extrapolation techniques. But not all implicit methods are suitable. A good method must have a stability region ${\displaystyle S}$ covering all or most of the negative half-plane ${\displaystyle {\mathbb {C} }^{-}}$.

The stability region is determined by applying the method to the Dahlquist test equation²¹ ${\displaystyle \,{\dot {u}}=\lambda u}$, and the stability region ${\displaystyle S}$ is the set of ${\displaystyle h\lambda \in {\mathbb {C} }}$ for which the method produces bounded solutions. A method with ${\displaystyle S\supset {\mathbb {C} }^{-}}$, i.e., where ${\displaystyle S}$ contains the entire left half-plane, is called A-stable.

There are many A-stable Runge-Kutta methods of high order to choose from, but unfortunately, for A-stable linear multistep methods the convergence order is limited to ${\displaystyle p=2}$. Therefore one usually has to settle for less. Among multistep methods the best choice for stiff equations are the backward differentiation (BDF) methods, to which the implicit Euler method belongs. There are BDF methods of convergence orders ${\displaystyle p=1:6}$, but only orders ${\displaystyle p\leq 2}$ are A-stable. While the stability region remains large for higher orders, it eventually deteriorates, and only methods up to order ${\displaystyle p=5}$ are regularly used.

Software based on BDF-type methods (or similar) are Matlab's ode15s, and C or Fortran codes such as CVODE, LSODE, MEBDF, DASSL. The software is elaborate, robust and of high complexity, uses variable order and variable time step adaptivity, and usually offers many options for problems with special structure, or for special handling of Jacobian matrices. Accuracy criteria can be specified both for absolute and relative error estimates. Some codes include options for nonstiff equations, or for differential-algebraic equations.

For implicit Runge-Kutta methods applied to the test equation, the differential equation is replaced by a recursion ${\displaystyle u_{n+1}=R(h\lambda )\,u_{n}}$, where the rational function ${\displaystyle R}$ is referred to as the stability function. If

${\displaystyle {\mathrm {Re} }\,z\leq 0\,\Rightarrow \,|R(z)|\leq 1\,,}$

the method is A-stable. There are A-stable methods of high order, but few implementations. An advantage of A-stable Runge-Kutta methods is that

${\displaystyle M_{2}[hA]\leq 0\,\Rightarrow \,\|R(hA)\|_{2}\leq 1\,,}$

i.e., if ${\displaystyle A}$ is negative definite, then the method produces a contractive recursion. (This is a variant of the von Neumann inequality.) Unfortunately, this only carries over to nonlinear problems under additional conditions. B-stable Runge-Kutta methods are a subset of A-stable methods, such that if ${\displaystyle \,M_{2}[hf]\leq 0\,}$ in a nonlinear problem, then the method produces a contractive recursion.²²²³²⁴ Other special stability requirements are also common, such as L-stability.²⁵ This is another subset of A-stable method, with the additional requirement ${\displaystyle \,R(\infty )=0\,}$ for improved damping of eigenmodes for which ${\displaystyle \,{\mathrm {Re} }\,h\lambda \ll -1}$.

Because implicit Runge-Kutta methods have high computational complexity, special efficiency requirements are also considered when software is developed. Among efficient Runge-Kutta software for stiff problems are Matlab's ode23s of orders ${\displaystyle 2}$ and ${\displaystyle 3}$, and the Fortran code RADAU5²⁶, implementing the B- and L-stable ${\displaystyle 5^{\mathrm {th} }}$ order Radau IIa method. For a general survey of diagonally implicit Runge-Kutta methods used for stiff equations, see Kennedy and Carpenter.²⁷

In all cases, the practical numerical solution of stiff equations needs dedicated, professional software. Even though it is a routine task, success may often require a careful study of the suitable settings of the codes and a proper understanding of the problem to be solved.

1. Is stiffness fully understood? In the literature it is sometimes suggested that no precise definition of stiffness exists. This appears to be due to the widespread mention or use of the stiffness ratio, whose obvious shortcomings are long since recognized. Sometimes operational criteria such as the actual discretization method and accuracy requirements are brought in to resolve the issues. But this over-complicates what is basically a simple issue. Thus, to the practitioner it is obvious that there is a distinction between nonstiff and stiff equations; it must therefore be possible to describe this distinction mathematically. The stiffness indicator provides a simple necessary criterion characterizing and quantifying this distinction. Today it is safe to say that stiffness is a complex but fully understood phenomenon, with excellent, efficient, and reliable dedicated software available.

2. Stiffness indicator. The stiffness indicator ${\displaystyle s_{2}[f']}$ is more robust than the scaled divergence ${\displaystyle \tau [f']}$. While the former only involves the two extreme logarithmic values, the latter makes use of the intermediate logarithmic values as well, in spite of the fact that they do not have a similar impact on stability. The advantage of the scaled divergence is that it bypasses eigenvalue computations and supports an a priori analytical understanding of the dynamics, as illustrated in the examples above. An approximate quantification is usually sufficient, as stiffness typically is a matter of orders of magnitude.

3. Dissipative vs. conservative systems. The notion of stiffness as discussed above is a measure of dissipativity, i.e., of damping or energy loss. Some sources suggest that there are also stiff conservative systems, in particular if derived from method-of-lines discretizations of hyperbolic PDEs. But such systems belong in a different category, of (potentially) highly oscillatory systems. With eigenvalues located far out on the imaginary axis, the corresponding high frequencies can be suppressed by choosing an implicit method using large steps. However, in accordance with the sampling theorem, high frequency phenomena are then no longer resolved and wave forms may be distorted. The distinction between dissipative and conservative problems is significant, and the computational approach is different. While stiffness is closely related to the irreversibility of parabolic problems, hyperbolic problems have no intrinsic damping and feature invariants, such as conserving energy. Similarly, separable Hamiltonian systems have a divergence free vector field and hence are "area preserving" (conserve phase volume). Due to this structure, the stiffness indicator is zero.

4. Etymology. The term "stiff" was introduced by Curtiss and Hirschfelder. According to Hirschfelder,²⁸ the term was chosen because the first examples were associated with servo systems, in which there was a "tight coupling" between the servo and the controlled system. A similar effect can also be seen in control systems with a high gain negative feedback, where higher gain increases stiffness. Another suggested connection is the (mechanical) notion of stiffness, as represented by the second order equation

${\displaystyle m{\ddot {x}}+c{\dot {x}}+kx=F\,,}$

where ${\displaystyle m}$ represents mass, ${\displaystyle c}$ the damping coefficient, ${\displaystyle k}$ the spring constant, and ${\displaystyle F}$ an external applied force. The idea is that equations with "stiff springs" (large ${\displaystyle k}$) cause stiffness in the mathematical sense. Taking ${\displaystyle m=1}$, the equation is rewritten as

${\displaystyle {\begin{aligned}{\dot {x}}&=y\\{\dot {y}}&=-kx-cy+F\,,\end{aligned}}}$

which is a system of the form ${\displaystyle {\dot {u}}=Au+g}$. The eigenvalues of ${\displaystyle A}$ are

${\displaystyle \lambda _{1,2}\,=\,{\frac {-c\pm {\sqrt {c^{2}-4k}}}{2}}\,;}$

these are claimed to be large if ${\displaystyle k}$ is large. For critical damping, the parameters must satisfy ${\displaystyle c=2{\sqrt {mk}}}$, which with ${\displaystyle m=1}$ leads to

${\displaystyle \lambda _{1,2}\,=\,-c/2\,.}$

Interestingly, computing the stiffness indicator (here identical to the scaled divergence) of the vector field one obtains

${\displaystyle s_{2}[A]\,=\,\tau [A]\,=\,-c/2\,,}$

which is large if and only if the damping constant ${\displaystyle c}$ is large, independent of the spring constant. Thus the single cause of stiffness is the damper, which dissipates energy. It is difficult to avoid the conclusion that the term "stiff" is a misnomer, which quickly got established in a new context with a different meaning. While in mechanical engineering a large spring constant typically goes with a large damping constant in order to have (near) critical damping, the confusion can be understood.

See also

See also

• Backward differentiation formula, a family of implicit methods especially used for the solution of stiff differential equations
• Condition number
• Differential inclusion, an extension of the notion of differential equation that allows discontinuities, in part as way to sidestep some stiffness issues
• Explicit and implicit methods

Notes

Notes

1. G Söderlind, LO Jay, M Calvo (2015). "Stiffness 1952--2012: Sixty years in search of a definition." BIT 55, pp 531-558.
2. G Söderlind (2024). "Logarithmic Norms." Berlin-Heidelberg-New York: Springer Series in Computational Mathematics SCM vol 63.
3. E. Hairer, G. Wanner (1996). "Solving orindary differential equations II. Stiff and differential-algebraic problems", 2nd ed. Berlin-Heidelberg-New York: Springer
4. C.F. Curtiss, J.O. Hirschfelder (1952). "Integration of stiff equations". Proc. Nat. Acad. Sci. 38, pp 235-243.
5. A. Prothero, A. Robinson (1974). "On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations". Math. Comp. 28, 145-162.
6. J.D. Lambert (1992). "Numerical Methods for Ordinary Differential Systems", pp 216-217. New York: Wiley, ISBN 978-0-471-92990-1.
7. G.D. Byrne and A.C. Hindmarsh (1987). "Stiff ODE Solvers: A review of current and coming attractions", p. 3ff. J. Comp. Phys. 70, 1-62.
8. J.D. Lambert (1973). "Cmmputational methods in ordinary differential equations", p. 232. London: John Wiley & Sons.
9. S. Artemiev, T. Averina (1997). "Numerical analysis of systems of ordinary and stochastic differential equations", p. 6. Utrecht: VSP.
10. K. Dekker, J.G. Verwer (1984). "Stability of Runge-Kutta methods for stiff nonlinear differential equations", p. 5. New York: North Holland.
11. L.F. Shampine (1985). "What is Stiffness?" In: R.C. Aiken (ed.), "Stiff Computation", p. 4. New York: Oxford University Press
12. G. Söderlind, L.O. Jay, M. Calvo (2015). "Stiffness 1952-2012. Sixty years in search of a definition. BIT Numerical Mathematics 55, 531-558.
13. ibid.
14. G. Söderlind (2024). "Logarithmic Norms". Springer Series in Computational Mathematics vol 63.
15. G. Söderlind, L.O. Jay, M. Calvo (2015). "Stiffness 1952-2012. Sixty years in search of a definition. BIT Numerical Mathematics 55, 531-558.
16. G. Söderlind (2024). "Logarithmic Norms", Ch. 15. Springer Series in Computational Mathematics, vol 63.
17. F. Mazzia, F. Iavernaro (2003). "Test Set for Initial Value Problem Solvers". Department of Mathematics, University of Bari, ITALY. http://archimede.dm.uniba.it/~testset/CWI_reports/testset2003r22.pdf
18. L.F. Shampine. https://es.mathworks.com/company/newsletters/articles/stiff-differential-equations.html
19. G Söderlind, LO Jay, M Calvo (2015). "Stiffness 1952--2012: Sixty years in search of a definition." BIT 55, pp 531-558.
20. G Söderlind (2024). Logarithmic Norms, Ch. 21. Springer Series in Computational Mathematics vol 63.
21. G. Dahlquist (1963). "A special stability problem for linear multistep methods". BIT 3, 27–43, doi:10.1007/BF01963532, hdl:10338.dmlcz/103497, S2CID 120241743
22. J.C. Butcher (1975). "A stability property of implicit Runge–Kutta methods". BIT 15, 358–361
23. J.C. Butcher (2008). "Numerical Methods for Ordinary Differential Equations", 2nd ed. New York: Wiley
24. E. Hairer, E., G. Wanner (1991). "Solving Ordinary Differential Equations II. Stiff and Differential– Algebraic Problems". Berlin–Heidelberg–New York: Springer.
25. Ehle (1969).
26. E. Hairer, E., G. Wanner (1991). "Solving Ordinary Differential Equations II. Stiff and Differential– Algebraic Problems". Berlin–Heidelberg–New York: Springer.
27. C.A. Kennedy, M.H. Carpenter (2016). "Diagonally Implicit Runge–Kutta Methods for Ordinary Differential Equations. A Review". NASA/TM-2016-219173, https://ntrs.nasa.gov/citations/20160005923
28. J.O. Hirshfelder (1963). "Applied Mathematics as used in Theoretical Chemistry". American Mathematical Society Symposium: 367–376.

References

References

• Burden, Richard L.; Faires, J. Douglas (1993), Numerical Analysis (5th ed.), Boston: Prindle, Weber and Schmidt, ISBN 0-534-93219-3.
• Dahlquist, Germund (1963), "A special stability problem for linear multistep methods", BIT, 3 (1): 27–43, doi:10.1007/BF01963532, hdl:10338.dmlcz/103497, S2CID 120241743.
• Eberly, David (2008), Stability analysis for systems of differential equations (PDF).
• Ehle, B. L. (1969), On Padé approximations to the exponential function and A-stable methods for the numerical solution of initial value problems (PDF), University of Waterloo.
• Gear, C. W. (1971), Numerical Initial-Value Problems in Ordinary Differential Equations, Englewood Cliffs: Prentice Hall, Bibcode:1971nivp.book.....G.
• Gear, C. W. (1981), "Numerical solution of ordinary differential equations: Is there anything left to do?", SIAM Review, 23 (1): 10–24, doi:10.1137/1023002.
• Hairer, Ernst; Wanner, Gerhard (1996), Solving ordinary differential equations II: Stiff and differential-algebraic problems (second ed.), Berlin: Springer-Verlag, ISBN 978-3-540-60452-5.
• Hirshfelder, J. O. (1963), "Applied Mathematics as used in Theoretical Chemistry", American Mathematical Society Symposium: 367–376.
• Iserles, Arieh; Nørsett, Syvert (1991), Order Stars, Chapman & Hall, ISBN 978-0-412-35260-7.
• Kreyszig, Erwin (1972), Advanced Engineering Mathematics (3rd ed.), New York: Wiley, ISBN 0-471-50728-8.
• Lambert, J. D. (1977), D. Jacobs (ed.), "The initial value problem for ordinary differential equations", The State of the Art in Numerical Analysis, New York: Academic Press: 451–501.
• Lambert, J. D. (1992), Numerical Methods for Ordinary Differential Systems, New York: Wiley, ISBN 978-0-471-92990-1.
• Mathews, John; Fink, Kurtis (1992), Numerical methods using MATLAB.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Section 17.5. Stiff Sets of Equations". Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8. Archived from the original on 2011-08-11. Retrieved 2011-08-17.
• Shampine, L. F.; Gear, C. W. (1979), "A user's view of solving stiff ordinary differential equations", SIAM Review, 21 (1): 1–17, doi:10.1137/1021001.
• Wanner, Gerhard; Hairer, Ernst; Nørsett, Syvert (1978), "Order stars and stability theory", BIT, 18 (4): 475–489, doi:10.1007/BF01932026, S2CID 8824105.
• Stability of Runge–Kutta Methods [1]

External links

External links

• An Introduction to Physically Based Modeling: Energy Functions and Stiffness
• Stiff systems Lawrence F. Shampine and Skip Thompson Scholarpedia, 2(3):2855. doi:10.4249/scholarpedia.2855