Stiffness Detection

Overview

Many differential equations exhibit some form of stiffness, which restricts the step size and hence effectiveness of explicit solution methods.

A number of implicit methods have been developed over the years to circumvent this problem.

For the same step size, implicit methods can be substantially less efficient than explicit methods, due to the overhead associated with the intrinsic linear algebra.

This cost can be offset by the fact that, in certain regions, implicit methods can take substantially larger step sizes.

Several attempts have been made to provide user-friendly codes that automatically attempt to detect stiffness at runtime and switch between appropriate methods as necessary.

A number of strategies that have been proposed to automatically equip a code with a stiffness detection device are outlined here.

Particular attention is given to the problem of estimation of the dominant eigenvalue of a matrix in order to describe how stiffness detection is implemented in (RefLink[NDSolve, paclet:ref/NDSolve]).

Numerical examples illustrate the effectiveness of the strategy.

Initialization

Load some packages with predefined examples and utility functions:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Needs["FunctionApproximations`"];

Introduction

Consider the numerical solution of initial value problems:

y'(t) == f(t, y(t)), y(0) == (Subscript[y, 0]), f: ℝ ×  (ℝ^n) ↦  (ℝ^n).

Stiffness is a combination of problem, solution method, initial condition and local error tolerances.

Stiffness limits the effectiveness of explicit solution methods due to restrictions on the size of steps that can be taken.

Stiffness arises in many practical systems as well as in the numerical solution of partial differential equations by the method of lines.

Example

The Van der Pol oscillator is a non-conservative oscillator with nonlinear damping and is an example of a stiff system of ordinary differential equations:

with ε==(3/1000);

consider initial conditions

and solve over the interval t∈[0,10].

The method "StiffnessSwitching" uses a pair of extrapolation methods by default:

Solution

This loads the problem from a package:

system = GetNDSolveProblem["VanderPol"];

Solve the system numerically using a nonstiff method:

solns = NDSolve[system, {T, 0, 10}, Method -> "Extrapolation"];

Solve the system using a method that switches when stiffness occurs:

sols = NDSolve[system, {T, 0, 10}, 
Method -> {"StiffnessSwitching","NonstiffTest"->False}];

Here is a plot of the two solution components. The sharp peaks (in blue) extend out to about 450 in magnitude and have been cropped:

Plot[Evaluate[Part[sols, 1, All, 2]], {T, 0, 10}, 
PlotStyle -> {{Red}, {Blue}}, Axes -> False, Frame -> True]

Stiffness can often occur in regions that follow rapid transients.

This plots the step sizes taken against time:

StepDataPlot[sols]

The problem is that when the solution is changing rapidly, there is little point using a stiff solver, since local accuracy is the dominant issue.

For efficiency, it would be useful if the method could automatically detect regions where local accuracy (and not stability) is important.

Linear Stability

Linear stability theory arises from the study of Dahlquist's scalar linear test equation:

y'(t) == λ  y(t),    λ ∈ ℂ,    Re(λ) < 0

as a simplified model for studying the initial value problem ().

Stability is characterized by analyzing a method applied to () to obtain

(Subscript[y, n+1])== R(z)  (Subscript[y, n])

where z==h λ and R (z) is the (rational) stability function.

The boundary of absolute stability is obtained by considering the region:

Explicit Euler Method

The explicit or forward Euler method

applied to () gives

The shaded region represents instability, where (Abs[R(z)])>1.

OrderStarPlot[1+z,1,z,FrameTicks->True]

The linear stability boundary (LSB) is often taken as the intersection with the negative real axis.

For the explicit Euler method LSB==–2.

For an eigenvalue of λ==-1, linear stability requirements mean that the step size needs to satisfy h<2, which is a very mild restriction.

However, for an eigenvalue of  λ==-(10^6), linear stability requirements mean that the step size needs to satisfy h < 2×(10^(-6)), which is a very severe restriction.

Example

This example shows the effect of stiffness on the step size sequence when using an explicit Runge-Kutta method to solve a stiff system.

This system models a chemical reaction:

system = GetNDSolveProblem["Robertson"];

The system is solved by disabling the built-in stiffness detection:

sol = NDSolve[system, Method->{"ExplicitRungeKutta", "StiffnessTest"->False}];

The step size sequence starts to oscillate when the stability boundary is reached:

StepDataPlot[sol]

The built-in detection does an excellent job of locating when stiffness occurs:

sol = NDSolve[system, Method->{"ExplicitRungeKutta", "StiffnessTest"->True}];

Implicit Euler Method

The implicit or backward Euler method:

applied to () gives:

The method is unconditionally stable for the entire left half-plane.

OrderStarPlot[1/(1-z),1,z,FrameTicks->True]

This means that to maintain stability there is no longer a restriction on the step size.

The drawback is that an implicit system of equations now has to be solved at each integration step.

Type Insensitivity

A type-insensitive solver recognizes and responds efficiently to stiffness at each step and so is insensitive to the (possibly changing) type of the problem.

One of the most established solvers of this class is LSODA [H83], [P83].

Later generations of LSODA such as CVODE no longer incorporate a stiffness detection device. The reason is because LSODA uses norm bounds to estimate the dominant eigenvalue and these bounds, as will be seen later, can be quite inaccurate.

The low order of A(α)-stable BDF methods means that LSODA and CVODE are not very suitable for solving systems with high accuracy or systems where the dominant eigenvalue has a large imaginary part. Alternative methods, such as those based on extrapolation of linearly implicit schemes, do not suffer from these issues.

Much of the work on stiffness detection was carried out in the 1980s and 1990s using standalone FORTRAN codes.

New linear algebra techniques and efficient software have since become available and these are readily accessible in the Wolfram Language.

Stiffness can be a transient phenomenon, so detecting nonstiffness is equally important [S77], [B90].

"StiffnessTest" Method Option

There are several approaches that can be used to switch from a nonstiff to a stiff solver.

Direct Estimation

A convenient way of detecting stiffness is to directly estimate the dominant eigenvalue of the Jacobian J of the problem (see [S77], [P83], [S83], [S84a], [S84c], [R87] and [HW96]).

Such an estimate is often available as a by-product of the numerical integration and so it is reasonably inexpensive.

If v denotes an approximation to the eigenvector corresponding to the dominant eigenvalue of the Jacobian, with ‖v‖ sufficiently small, then by the mean value theorem a good approximation to the leading eigenvalue is

Richardson's extrapolation provides a sequence of refinements that yields a quantity of this form, as do certain explicit Runge–Kutta methods.

Cost is at most two function evaluations, but often at least one of these is available as a by-product of the numerical integration, so it is reasonably inexpensive.

Let LSB denote the linear stability boundary—the intersection of the linear stability region with the negative real axis.

The product h(Overscript[λ, ∼]) gives an estimate that can be compared to the linear stability boundary of a method in order to detect stiffness:

|h(Overscript[λ, ∼])|≤ s|LSB|

where s is a safety factor.

Description

The methods "DoubleStep", "Extrapolation", and "ExplicitRungeKutta" have the option "StiffnessTest", which can be used to identify whether the method applied with the specified (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal]) and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal]) tolerances to a given problem is stiff.

The method option "StiffnessTest" itself accepts a number of options that implement a weak form of () where the test is allowed to fail a specified number of times.

The reason for this is that some problems can be only mildly stiff in a certain region and an explicit integration method may still be efficient.

"NonstiffTest" Method Option

The "StiffnessSwitching" method has the option "NonstiffTest", which is used to switch back from a stiff method to a nonstiff method.

The following settings are allowed for the option "NonstiffTest":

Switching to a Nonstiff Solver

An approach that is independent of the stiff method is used.

Given the Jacobian J (or an approximation), compute one of:

Norm Bound: ‖ J ‖

Spectral Radius: ρ(J)==max (Abs[(Subscript[λ, i])])

Dominant Eigenvalue (Subscript[λ, i]) : |(Subscript[λ, i])| > | (Subscript[λ, j]) |.

Many linear algebra techniques focus on solving a single problem to high accuracy.

For stiffness detection, a succession of problems with solutions to one or two digits is adequate.

For a numerical discretization

consider a sequence k of matrices in some sub-interval(s)

The spectra of the succession of matrices often change very slowly from step to step.

The goal is to find a way of estimating (bounds on) dominant eigenvalues of a succession of matrices (Subscript[J, (Subscript[t, i])])that:

NormBound

A simple and efficient technique of obtaining a bound on the dominant eigenvalue is to use the norm of the Jacobian (Norm2[J, p]) where typically p == 1 or p == ∞.

The method has complexity O((n^2)), which is less than the work carried out in the stiff solver.

This is the approach used by LSODA.

The setting "NormBound" of the option "NonstiffTest" computes (Norm2[J, 1]) and (Norm2[J, ∞]) and returns the smaller of the two values.

Example

The following Jacobian matrix arises in the numerical solution of the Van der Pol system using a stiff solver:

a = {{0.,1.},{2623.532160943381,-69.56342161343568}};

Bounds based on norms overestimate the spectral radius by more than an order of magnitude:

{Abs[First[Eigenvalues[a]]], Norm[a, 1], Norm[a, Infinity]}

Direct Eigenvalue Computation

For small problems (n≤32) it can be efficient just to compute the dominant eigenvalue directly.

The setting "Direct" of the option "NonstiffTest" computes the dominant eigenvalue of J using the same LAPACK routines as Eigenvalues.

For larger problems the cost of direct eigenvalue computation is O((n^3)), which becomes prohibitive when compared to the cost of the linear algebra work in a stiff solver.

A number of iterative schemes have been implemented for this purpose. These effectively work by approximating the dominant eigenspace in a smaller subspace and using dense eigenvalue methods for the smaller problem.

The Power Method

Shampine has proposed the use of the power method for estimating the dominant eigenvalue of the Jacobian [S91].

The power method is perhaps not a very well-respected method, but has received a resurgence of interest due to its use in Google's page ranking.

The power method can be used when

Description

Given a starting vector (Subscript[v, 0])∈(ℂ^n), compute

The Rayleigh quotient is used to compute an approximation to the dominant eigenvalue:

In practice, the approximate eigenvector is scaled at each step:

Properties

The power method converges linearly with rate

which can be slow.

In particular, the method does not converge when applied to a matrix with a dominant complex conjugate pair of eigenvalues.

Generalizations

The power method can be adapted to overcome the issue of equimodular eigenvalues (e.g. NAPACK).

However the modification does not generally address the issue of the slow rate of convergence for clustered eigenvalues.

There are two main approaches to generalizing the power method:

Although the methods work quite differently, there are a number of core components that can be shared and optimized.

Subspace and Krylov iteration cost O(m (n^2)) operations.

They project an n×n matrix to an m×m matrix, where m<<n .

The small matrix represents the dominant eigenspace and approximation uses dense eigenvalue routines.

Subspace Iteration

Subspace (or simultaneous) iteration generalizes the ideas in the power method by acting on m vectors at each step.

Start with an orthonormal set of vectors (V^((0))) ==(ℂ^(n × m)), where usually m<<n :

Form the product with the matrix A:

In order to prevent all vectors from converging to multiples of the same dominant eigenvector (Subscript[v, 1]) of A, they are orthonormalized:

The orthonormalization step is expensive compared to the matrix product.

Rayleigh–Ritz Projection

Input: matrix A and an orthonormal set of vectors V

The matrix S has small dimension m×m.

Note that the Schur decomposition can be computed in real arithmetic when S∈(ℝ^(m×m)) using a quasi upper triangular matrix T.

Convergence

Subspace (or simultaneous) iteration generalizes the ideas in the power method by acting on m vectors at each step.

SRRIT converges linearly with rate

In particular, the rate for the dominant eigenvalue is

Therefore it can be beneficial to take, for example, m==3 or more even if we are only interested in the dominant eigenvalue.

Error Control

A relative error test on successive approximations' dominant eigenvalue is

This is not sufficient since it can be satisfied when convergence is slow.

If |(Subscript[λ, i])| ==|(Subscript[λ, i-1])| or |(Subscript[λ, i])| ==|(Subscript[λ, i+1])| then the i(^th) column of (Q^((k))) is not uniquely determined.

The residual test used in SRRIT is:

where ((OverHat[Q])^((k)))==(Q^((k))) (U^((k))), (Subscript[(OverHat[q]), i]^((k))) is the i(^th) column of ((OverHat[Q])^((k))), and (Subscript[t, i]^((k))) is the i(^th) column of (T^((k))).

This is advantageous since it works for equimodular eigenvalues.

The first column position of the upper triangular matrix (T^((k))) is tested because of the use of an ordered Schur decomposition.

Implementation

There are several implementations of subspace iteration.

The implementation for use in "NonstiffTest" is based on:

"An attractive feature of SRRIT is that it displays monotonic consistency, that is, as the convergence tolerance decreases so does the size of the computed residuals" [LS96].

SRRIT makes use of an ordered Schur decomposition where eigenvalues of the largest modulus appear in the upper-left entries.

Modified Gram–Schmidt with reorthonormalization is used to form (Q^((k))), which is faster than Householder transformations.

The approximate dominant subspace (Subscript[V, (Subscript[t, i])]^((k))) at integration time (Subscript[t, i]) is used to start the iteration at the next integration step (Subscript[t, i+1]):

KrylovIteration

Given an n×m matrix V whose columns (Subscript[v, i]) comprise an orthogonal basis of a given subspace 𝒱:

The Rayleigh–Ritz procedure consists of computing H==(V^T)A V and solving the associated eigenproblem H (Subscript[y, i])==(Subscript[θ, i])(Subscript[y, i]).

The approximate eigenpairs of the original problem (Subscript[(OverTilde[λ]), i]), (Subscript[(OverTilde[x]), i]) satisfy (OverTilde[λ])==(Subscript[θ, i]) and (Subscript[(OverTilde[x]), i]) == V (Subscript[y, i]), which are called Ritz values and Ritz vectors.

The process works best when the subspace 𝒱 approximates an invariant subspace of A.

This process is effective when 𝒱 is equal to the Krylov subspace associated with a matrix A and a given initial vector x as

Description

The method of Arnoldi is a Krylov-based projection algorithm that computes an orthogonal basis of the Krylov subspace and produces a projected m×m matrix H with m<<n.

Input: matrix A, the number of steps m, an initial vector (Subscript[v, 1]) of norm 1.

Output: ((Subscript[V, m]),(Subscript[H, m]),f,β) with β==(Subscript[‖f‖, 2])

In the case of Arnoldi, H has an unreduced upper Hessenberg form (upper triangular with an additional nonzero subdiagonal).

Orthogonalization is usually carried out by means of a Gram–Schmidt procedure.

The quantities computed by the algorithm satisfy

The residual f gives an indication of proximity to an invariant subspace and the associated norm β indicates the accuracy of the computed Ritz pairs:

Restarting

The Ritz pairs converge quickly if the initial vector x is rich in the direction of the desired eigenvalues.

When this is not the case then a restarting strategy is required in order to avoid excessive growth in both work and memory.

There are several strategies for restarting, in particular:

Explicit restart is relatively simple to implement, but implicit restart is more efficient since it retains the relevant eigeninformation of the larger problem. However implicit restart is difficult to implement in a numerically stable way.

An alternative that is much simpler to implement, but achieves the same effect as implicit restart, is a Krylov–Schur method [S01].

Implementation

A number of software implementations are available, in particular:

The implementation in "NonstiffTest" is based on Krylov–Schur iteration.

Automatic Strategy

The "Automatic" setting uses an amalgamation of the methods as follows.

Step Rejections

Caching of the time of evaluation ensures that the dominant eigenvalue estimate is not recomputed for rejected steps.

Stiffness detection is also performed for rejected steps since:

Iterative Method Options

The iterative methods of "NonstiffTest" have options that can be modified:

Options[NDSolve`SubspaceIteration]

Options[NDSolve`KrylovIteration]

The default tolerance aims for just one correct digit, but often obtains substantially more accurate values—especially after a few successful iterations at successive steps.

The default values limiting the number of iterations are:

If these values are set too large then a convergence failure becomes too costly.

In difficult problems, it is better to share the work of convergence across steps. Since the methods effectively refine the basis vectors from the previous step, there is a reasonable chance of convergence in subsequent steps.

Latency and Switching

It is important to incorporate some form of latency in order to avoid a cycle where the "StiffnessSwitching" method continually tries to switch between stiff and nonstiff methods.

The options "MaxRepetitions" and "SafetyFactor" of "StiffnessTest" and "NonstiffTest" are used for this purpose.

The default settings allow switching to be quite reactive, which is appropriate for one-step integration methods.

Examples

Van der Pol

Select an example system:

system = GetNDSolveProblem["VanderPol"];

StiffnessTest

The system is integrated successfully with the given method and the default option settings for "StiffnessTest":

NDSolve[system, Method -> "ExplicitRungeKutta"]

A longer integration is aborted and a message is issued when the stiffness test condition is not satisfied:

NDSolve[system, {T, 0, 10}, Method -> "ExplicitRungeKutta"]

Using a unit safety factor and specifying that only one stiffness failure is allowed effectively gives a strict test. The specification uses the nested method option syntax:

NDSolve[system,Method -> {"ExplicitRungeKutta", "StiffnessTest"->{True, "MaxRepetitions" -> {1,1}, "SafetyFactor"->1} }]

NonstiffTest

For such a small system, direct eigenvalue computation is used.

The example serves as a good test that the overall stiffness switching framework is behaving as expected.

Set up a function to monitor the switch between stiff and nonstiff methods and the step size taken. Data for the stiff and nonstiff solvers is put in separate lists by using a different tag for (RefLink[Sow, paclet:ref/Sow]):

SetAttributes[SowSwitchingData, HoldFirst];

SowSwitchingData[told_,t_,method_NDSolve`StiffnessSwitching]:=
(Sow[{t,t-told},method["ActiveMethodPosition"]];
told=t;);

Solve the system and collect the data for the method switching:

T0 = 0;
data=
Last[
Reap[
sol=NDSolve[system,{T,0,10},
Method->StiffnessSwitching,
"MethodMonitor":>(SowSwitchingData[T0, T,NDSolve`Self];)
];
]
];

Plot the step sizes taken using an explicit solver (blue) and an implicit solver (red):

ListLogPlot[data,Axes->False,Frame->True, PlotStyle->{Blue,Red}]

Compute the number of nonstiff and stiff steps taken (including rejected steps):

Map[Length,data]

CUSP

The cusp catastrophe model for the nerve impulse mechanism [Z72]:

Combining with the Van der Pol oscillator gives rise to the CUSP system [HW96]:

where

and σ==(1/144) and ε==(1/(10^4)).

Discretization of the diffusion terms using the method of lines is used to obtain a system of ODEs of dimension 3 n==96.

Unlike the Van der Pol system, because of the size of the problem, iterative methods are used for eigenvalue estimation.

Step Size and Order Selection

Select the problem to solve:

system = GetNDSolveProblem["CUSP-Discretized"];

Set up a function to monitor the type of method used and step size. Additionally, the order of the method is included as a tooltip:

SetAttributes[SowOrderData, HoldFirst];

SowOrderData[told_,t_,method_NDSolve`StiffnessSwitching]:=
(Sow[
Tooltip[{t,t-told},method["DifferenceOrder"]],
method["ActiveMethodPosition"]
];
told=t;);

Collect the data for the order of the method as the integration proceeds:

T0 = 0;
data=
Last[
Reap[
sol=NDSolve[system,
Method->"StiffnessSwitching",
"MethodMonitor":>(SowOrderData[T0, T, NDSolve`Self];)
];
]
];

Plot the step sizes taken using an explicit solver (blue) and an implicit solver (red). A tooltip shows the order of the method at each step:

ListLogPlot[data,Axes->False,Frame->True, PlotStyle->{Blue,Red}]

Compute the total number of nonstiff and stiff steps taken (including rejected steps):

Map[Length,data]

Jacobian Example

Define a function to collect the first few Jacobian matrices:

SetAttributes[StiffnessJacobianMonitor,HoldFirst];

StiffnessJacobianMonitor[i_,method_NDSolve`StiffnessSwitching]:=If[SameQ[method["ActiveMethodPosition"],2]&&i<5,
If[MatrixQ[#],
Sow[#];
i=i+1
]& @ method["Jacobian"]
];

i=0;
jacdata=Reap[sol=NDSolve[system,Method->"StiffnessSwitching","MethodMonitor":>(StiffnessJacobianMonitor[i,NDSolve`Self];)];
][[-1,1]];

A switch to a stiff method occurs near 0.00113425 and the first test for nonstiffness occurs at the next step (Subscript[t, k])≈0.00127887`.

Graphical illustration of the Jacobian (Subscript[J, (Subscript[t, k])]):

MatrixPlot[First[jacdata]]

Define a function to compute and display the first few eigenvalues of (Subscript[J, (Subscript[t, k])]), (Subscript[J, (Subscript[t, k+1])]), … and the norm bounds:

DisplayJacobianData[jdata_]:=
Module[{evdata, hlabels, vlabels},
evdata =
Map[Join[Eigenvalues[Normal[#], 4],{Norm[#,1],Norm[#,Infinity]}]&,jdata];
vlabels = {{""},{"(Subscript[λ, 1])"},{"(Subscript[λ, 2])"},{"(Subscript[λ, 3])"},{"(Subscript[λ, 4])"},{"\[LeftDoubleBracketingBar](Subscript[J, (Subscript[t, k])])(Subscript[\[RightDoubleBracketingBar], 1])"},{"\[LeftDoubleBracketingBar](Subscript[J, (Subscript[t, k])])(Subscript[\[RightDoubleBracketingBar], ∞])"}};
hlabels = Table[(Subscript[J, (Subscript[t, k])]),{k,Length[jdata]}];
Grid[MapThread[Join, {vlabels,Join[{hlabels},Transpose[evdata]]}],Frame->All]
];

DisplayJacobianData[jacdata]

Norm bounds are quite sharp in this example.

Korteweg–deVries

The Korteweg–deVries partial differential equation is a mathematical model of waves on shallow water surfaces:

We consider boundary conditions

and solve over the interval t∈[0,1].

Discretization using the method of lines is used to form a system of 192 ODEs.

Step Sizes

Select the problem to solve:

system = GetNDSolveProblem["Korteweg-deVries-PDE"];

The backward differentiation formula methods used in LSODA run into difficulties solving this problem:

First[Timing[sollsoda=NDSolve[system,Method->LSODA];]]

A plot shows that the step sizes rapidly decrease:

StepDataPlot[sollsoda]

In contrast, StiffnessSwitching immediately switches to using the linearly implicit Euler method, which needs very few integration steps:

First[Timing[sol = NDSolve[system,Method->"StiffnessSwitching"];]]

StepDataPlot[sol]

The extrapolation methods never switch back to a nonstiff solver once the stiff solver is chosen at the beginning of the integration.

Therefore this is a form of worst case example for the nonstiff detection.

Despite this, the cost of using subspace iteration is only a few percent of the total integration time.

Compute the time taken with switching to a nonstiff method disabled:

First[Timing[sol = NDSolve[system,
Method->{"StiffnessSwitching", "NonstiffTest"->False}];]]

Jacobian Example

Collect data for the first few Jacobian matrices using the previously defined monitor function:

i=0;
jacdata=Reap[sol=NDSolve[system,Method->"StiffnessSwitching","MethodMonitor":>(StiffnessJacobianMonitor[i,NDSolve`Self];)];
][[-1,1]];

Graphical illustration of the initial Jacobian (Subscript[J, (Subscript[t, 0])]):

MatrixPlot[First[jacdata]]

Compute and display the first few eigenvalues of (Subscript[J, (Subscript[t, k])]), (Subscript[J, (Subscript[t, k+1])]), … and the norm bounds:

DisplayJacobianData[jacdata]

Norm bounds overestimate slightly, but more importantly they give no indication of the relative size of real and imaginary parts.

Option Summary

StiffnessTest

Options of the method option "StiffnessTest".

NonstiffTest

Options of the method option "NonstiffTest".