Numerical instability and dynamical systems

Abstract

In philosophical studies regarding mathematical models of dynamical systems, instability due to sensitive dependence on initial conditions, on the one side, and instability due to sensitive dependence on model structure, on the other, have by now been extensively discussed. Yet there is a third kind of instability, which by contrast has thus far been rather overlooked, that is also a challenge for model predictions about dynamical systems. This is the numerical instability due to the employment of numerical methods involving a discretization process, where discretization is required to solve the differential equations of dynamical systems on a computer. We argue that the criteria for numerical stability, as usually provided by numerical analysis textbooks, are insufficient, and, after mentioning the promising development of backward analysis, we discuss to what extent, in practice, numerical instability can be controlled or avoided.

Appendix: A

This Appendix provides examples and mathematical details on some notions discussed in the paper, viz. stiff problems, zero-stability and absolute stability.

1.1 A.1 Stiff problems

Here is an one-dimensional example of stiff problem provided by Hairer and Wanner (1991):¹⁵

$$ {\frac{dy}{dx}=-50(y-\cos{x})} $$

(7)

On the left of Fig. 5 are the solution curves of the equation. They show a smooth solution, close to $y=\cos \limits {x}$, as well as a bunch of solutions that reach this smooth solution after a transient phase. However, some other numerical solutions do not behave the same way. On the right of Fig. 5 are two solutions with initial value y(0) = 1, and respectively step sizes h = 1.974/50 and 1.875/50. These solutions have no transient phase, and enable us to infer that, as soon as the step size is larger than 2/50, the numerical solution goes too far beyond the equilibrium and violent oscillations occur (Hairer and Wanner 1991, 2).

Fig. 5

On the left, solution curves with implicit Euler solution. On the right, explicit Euler solution for y(0) = 0, h = 1.974/50 and 1.875/50. (Hairer and Wanner 1991, 2)

Specific numerical methods have been developed in order to solve stiff problems (see Appendix Appendix).

1.2 A.2 Zero-stability

In order to understand the importance of a numerical method being zero-stable, let us take the example of a dynamical system defined by the following partial differential equation:

$$ {\frac{\partial u}{\partial t}=\frac{\partial u}{\partial x}} $$

(8)

with $u(x,0)=\cos \limits ^{2}x$ if $\lvert x \rvert \leq \pi /2$, u(x, 0) = 0 otherwise.

Let us take the leap frog numerical method, viz. the discretization $u^{n+1}_{j}=u^{n-1}_{j}+\lambda (u^{n}_{j+1}-u^{n}_{j-1})$, with the space step h = 0.04π and time step k = λh, where λ is the mesh ratio (see Trefethen 1994, 148). Figure 6 shows the two numerical solutions with respectively λ = 0.9 and λ = 1.1. As we can see, for λ < 1 there is a propagation of the wave as predicted by the theory. This is the regime where the numerical method is zero-stable. On the contrary, for λ > 1 numerical artefacts appear in the propagation of the wave. This is the regime where the numerical method is not zero-stable.

Fig. 6

Stable and unstable leap frog approximations to $\frac {\partial u}{\partial t}=\frac {\partial u}{\partial x}$. (Trefethen 1994, 149, Fig. 4.1.1)

As emphasized by Trefethen (1994), in the case where λ > 1: “The errors introduced at each step are not much bigger than before [i.e. for λ < 1] but they grow exponentially in subsequent time steps [...]. This rapid blow up of a sawtooth mode is typical of unstable finite difference formulas” (1994, chap 4, 148).

It is important to emphasize that the instability is not generated by round-off errors; it is generated by discretization errors which accumulate in the regime where the numerical method is not zero-stable. In both cases (λ < 1 and λ > 1) the numerical method is consistent. Nevertheless, the numerical solution to Eq. 8 would converge to the exact solution when h and k tend to zero in the regime of zero-stability, i.e. when λ < 1.¹⁶ Thus, as we have just seen, in order to avoid numerical instability, it can be crucial to investigate whether a numerical method is zero-stable.

Zero-stability can be applied to numerical methods used for ordinary differential equations. In those cases, it concerns the behavior of the numerical method for a fixed time value t, as time step h tends to zero (see Hairer et al. 1987, 378, and Fasshauer 2007, chap. 4). This is a particularly important condition for multi-step methods (Süli and Mayers; Hairer et al. 1987, 398; Trefethen 1994, chap. 1, 40 2003; Hairer et al. 1987, 398; Trefethen 1994, chap. 1, 40, 332). Zero-stability also concerns the numerical methods used for partial differential equations. In this case, it also concerns the behavior of the numerical method for a fixed time value t, as time step k tends to zero, but also as spatial step h tends to zero.

Consistency is a property required for a numerical method. A numerical method is said to be consistent if the local errors introduced at each step tend to zero when the discretization steps tend to zero. However, consistency is a necessary but not sufficient condition for solution convergence, which means that the discrete solution tends to the exact continuous one when the time step tends to zero. Zero-stability is also required: if a numerical method is consistent and zero-stable, then it is convergent. In other words, being zero-stable allows a numerical method to be convergent.¹⁷

1.3 A.3 Absolute stability and stiffness

In order to better understand what absolute stability is, and how it is connected with stiff problems, let us take the example of the initial value problem:

$$ \frac{dx}{dt}= -\mu x \ \ \ x(0)=a $$

(9)

The solution is an exponential function such as x(t) = ae^−μt with μ > 0.

The resolution of the initial value problem with the Euler forward method leads to the following discrete equation:

$$ x_{k+1}=x_{k}-h\mu x_{k}, \ \ \ x_{0}=a $$

(10)

The exact discrete solution is: x_k = a(1 − μh)^k. If $h<\frac {2}{\mu }$, both continuous and discrete solutions tend to zero when t or k tends to infinity. Nevertheless, when $h>\frac {2}{\mu }$, the discrete solution increases with k. This illustrates that the Euler method is not absolutely stable for all values of h. It is a conditionally stable numerical method.

Let us now use the Euler backward method in order to solve this same problem. In that case, the discrete equation becomes:

$$ x_{k+1}=x_{k}-h\mu x_{k+1}, \ \ \ x_{0}=a $$

(11)

The exact discrete solution is: x_k = x₀/(1 + μh)^k. Here, both continuous and discrete equations tend to zero when t and k tend to infinity, whatever the value of h is. This illustrates that the Euler backward method is an unconditionally stable numerical method.

The Euler forward method is an explicit method: one just needs to know the value x_k to directly infer the value x_k+ 1 at the next time step by applying the discrete time evolution. On the contrary, the Euler backward method is implicit: one cannot directly derive the value of x_k+ 1 with x_k, instead one generally has to solve a non-linear equation with x_k+ 1. Generally, explicit methods are not suitable to solve stiff problems – something which is sometimes part of the definition of stiffness itself! Thus, Hairer et al. (1987) assert that “stiff equations are problems for which explicit methods don’t work” (1987, 2). Implicit methods, on the contrary, are indeed generally suitable.¹⁸