EPODE

Stiff problems

The essence of the stiff phenomenon consists in the fact that the exact solution includes some components with a very fast decreasement that can be very hard to be followed by the numerical solution given by a step by step iterative process.

We mention some pragmatic definitions:

The stiff equations are the equations for which some implicit methods work better than the explicit ones.
A problem is stiff in a given integration interval if, for a given numerical code, the stepsize must be very strongly reduced.
The stiff differential equations are wrong conditioned in the computational sense.
An ordinary differential equation system is stiff in a given integration interval [0,T] if, in the exact solution, there is at least one component with a very large variation relative to the value 1/T.

Example 1. Let the scalar equation

y'(t)=f(t,y)=l y(t)+F'(t)-l F(t), t>= 0, l< < 0. The exact solution of the associated IV problem is

y(t)=F(t)+e^{{l t}[y₀-F(0)], We notice that, after a short time, the transitory component e^{{l t}[y₀-F(0)], also referred to as the stiff component, has no more a signifiant influence on the solution value. There is a subinterval in which the slow varying component, F(t), also referred to as the untransitory component or the smooth component, is prevalent in the exact solution. This subinterval covers almost all the integration interval. Applying the Euler's explicit rule with stepsize h to this problem, we get

y_n+1=y_n+hf(t_n,y_n)=(1+hl)[y_n-F(t_n)]+F(t_n)+hF'(t_n). If the numerical model is correct, the difference y_n-F(t_n) must decrease: this is true iff -2< hl< 0. Thus, the stepsize is retricted by some numerical stability condition.

In the transitory stage a problem can not be classified to be stiff one, only when the stiff component is negligible. An equation can not be referred to as stiff one (only the IVP associated to it), since the term is relative to the initial value, the integration interval and the error tolerance.

Example 2. Some components with different dials of variation can be meet especially in the case of the linear systems of the following form

y'(t)=Ay(t)+F(t), where A is a constant matrix with different dials values of the eigenvalues l_i. The above problem is stiff if

$ i: Re (l_i)< < 0;
$ j: |l_j|< < |l_i|;
not $ k: Re (l_k)> > 0;
not $ m: | Im (l_m)|> > 0, only if Re (l_m)< < 0.

We define the stiff ratio of the above system by

S(A)= max_{j=1, ... ,N}| Re (l_j)|/ min_{j=1, ... ,N}| Re (l_j)|.

The above system is stiff if Re (l_j)< 0, j=1, ... ,N, and S(A)> > 1.

Example 3. An important class of stiff problems with very large dimensions is derived from some partial differential equations. For example, let the problem

u_t=u_xx, t>0, 0<= x<= 1, u(t,0)=u(t,1)=0, t>0, u(0,x)=j(x) and u_i(t) the approximate values of the exact values u(t,x_i), produced by replacing u_xx with

d²u(t,x_i) / h²=[u(t,x_i+h)-2u(t,x_i)+u(t,x_i-h)]/ h², x_i=x_i-1+h, where h is the stepsize on x axis, h=1/m, m in N. Then u_i(t),i=1, ... ,m-1 are the components of the exact solution of the IV linear problem

y'(t)=Ay(t), t>0, A=m²tridiag(1,-2,1) The eigenvalues l_i of the matrix A,

l_i=m²[-2+2 cos(ip/ m)], i=1, ... ,m-1, are uniformly distributed in the interval (-4m²,0). When m increases (the division is more dense) the stiff ratio also increases.

Example 4. The singular perturbation problems are particular cases of stiff problems. The generic form of a such problem is

x'=f(t,x,y,e), e y'=g(t,x,y,e,) x(0)=e, y(0)=h, g(t,x,y,0)=0. In fact this system is a stiff one. For example, let f=g=x+y. Then the eigenvalues of the system are e^-1+O(1) and -1+O(e). When e decreases, the stiff ratio is increasing.

Such systems can be solved by some pseudo-stationary methods. We consider the reduced problem

x'=f(t,x,y,0), 0=g(t,x,y,0), x(0)=e. Unfortunately, there are some systems which can not be simplified in this manner since the presence of the control value e determines the evolution of the exact solution.

Example 5. If the system presents a variation on t, like

y'(t)=A(t)y(t), the affiliation to the stiff class is not easy to be establish. The problem is only local a stiff one. The variation of the eigenvalues of the Jacobian's matrix can not inform us about the perturbation propagation in the exact solution. In practice, a system is classified to be stiff one if, at least at a moment (at one t from the integration interval), the associated initial value problem is a locally stiff problem.

In physics and chemestry, any real system which is modelated by an ordinary differential equation system with components supposing some evolutions in different dials of time, also referred to as time constants, is treated as a stiff system. A large time constant determines the stiff character. The components of practical interest are those with a slow variation. Integrating the problem with an explicit method, the stepsize must be of small constant order.

Example 6. A nonlinear system is stiff if at least for one point t_s in [a,b] the associated linear system y'(t)=(d f/d y)(t_s,y(t_s))y(t), t>= t_s is stiff.

The meaning of the sign > > in the above mentionated inequalities is relative to the integration interval length and the error tolerance. An initial value problem can not be classified as stiff one, also when the eigenvalues are of very different values, if the integration interval can be compared with the transitory stage, or if the error tolerance e is not very small.

The stiff ratio of the above system is

S(A,a,b,e)=max_{i=1, ... ,N}{|l_i|(b-a)lg (1/ e),

where a,b are the extreme values of the integration interval and e is the error tolerance.

Depending on the stiff ratio the problems can be classified like in the following table:

Class	Stiff ratio
Stiff at limit	`S=O(10)`
Midle stiffly	`S=O(10²)`
Strongly stiffly	`O(10²)<= S <= O(10⁵)`
Extreme stiffly	`O(10⁶)<= S <= O(10⁸)`
Pathological stiffly	`O(10⁹)<=S`

In chemestry dynamics we can often meet a stiff ratio of order O(10⁶).

An ordinary differential equation system y'=f(t,y), numerical integrated on the interval [a,b], starting from the initial value y(a)=y_a, is stiff if the function f has a Lipschitz's constant very large relative to the interval length, and the error tolerance.

The principal characteristics of the stiff problems are:

the exact solutions are stable in the sense that small perturbations in the initial values are followed only by small perturbations in the exact solutions;
trying to solve the problem with the standard methods we get some strict restrictions on the stepsize from the stability conditions.

A problem necessarly must have these two characteristics to be declared a stiff one.