IF the step length used is too small, excessive computation time and round-off error result. We should also consider the opposite case, and ask whether there is any upper bound on step length. Often there is such a bound, and it is reached when the method becomes numerically unstable: the numerical solution produced no longer corresponds qualitatively with the exact solution because some bifurcation has occurred.
The traditional criterion for ensuring that a numerical method is stable is called absolute stability. Absolute-stability analysis of Runge--Kutta and other numerical methods is carried out using the linear model problem
where is complex. This has the analytical solution
The problem has a stable fixed point at y=0 for .
For systems of equations we generalize to the problem
where A is a matrix with distinct eigenvalues all lying in the negative half-plane so that again we have a stable fixed point at y=0.
We are interested in these linear model problems since they underlie the theory of the classification of fixed points. Since its eigenvalues are distinct, A is diagonalizable and we can reduce Eq.(43) to a set of independent equations of the form of Eq.(41). are then the eigenvalues of the Jacobian matrix. This is why we allow to be complex above; higher-dimensional systems may have complex eigenvalues, so we need to know the behaviour of the numerical solutions in complex -space.
The region of absolute stability for a method is then the set of values of h (real and non-negative) and (complex) for which as , i.e., for which the fixed point at the origin is stable. Thus we want the set of values of h and for which where S, the stability function, is the eigenvalue of the Jacobian of the Runge--Kutta map evaluated at the fixed point. In the case of systems, the modulus of the stability function is given by the spectral radius of the Jacobian of the Runge--Kutta map (the absolute value of the largest eigenvalue of the Jacobian) at the fixed point.
For example, let us derive the equation of the region of absolute stability for the family of Runge--Kutta methods derived earlier and presented in Eq.(25). Using the linear model problem we have
and the stability function is
Figure 1: The absolute stability regions of explicit p-stage, pth-order Runge--Kutta methods for are plotted in complex -space. The absolute stability regions are shown in grey. The ordinate and abscissa are and respectively. Notice that the size of the regions increases with the order of the method.
In general, recall that for an explicit method of order p, we wished to have Eq.(8) matching Eq.(4) up to terms of . Substituting and into Eq.(5):
Thus for an explicit p-stage method of order p (which is only possible for ), the stability function is
For an explicit q-stage method of order p<q,
where the s are determined by the particular Runge--Kutta method. In this case the free parameters may be used to maximize the area of the absolute-stability region. In practice, however, one does not get large increases in the area by doing this. It is interesting that in the optimal case p=q, the free parameters do not come into the expression for S, so all optimal methods of a given order will have the same absolute-stability region. We show the stability boundaries given by |S|=1 for in Fig. 1. Notice that the size of the absolute-stability region increases with the order of the method. Note also that h and will always be found as the pair in this model problem, due to the method of construction of S. It is clear that for a q-stage method, S is a polynomial of degree q in h. Since as , explicit methods all have bounded absolute-stability regions. Implicit Runge--Kutta methods have absolute-stability regions which can be unbounded. This is a great advantage in some situations, as we shall see below.