Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations

The front matter includes Table of Contents, List of Figures, List of Tables and Preface

Ordinary differential equations (ODEs) arise in many instances when using mathematical modeling techniques for describing phenomena in science, engineering, economics, etc. In most cases the model is too complex to allow one to find an exact solution or even an approximate solution by hand: an efficient, reliable computer simulation is required.

The term stability has been used in the literature for a large variety of different concepts. The basic, qualitative idea is that a model that produces a solution (output) for given data (input) should possess the property that if the input is perturbed by a small amount then the output should also be perturbed by only a small amount. But the precise application of this idea to initial value ordinary differential equations (ODEs), to boundary value ODEs, and to numerical methods has given rise to a multitude of definitions. The reader should therefore be careful, when speaking of stability, to distinguish between stability of problems and of numerical methods, and between stability of initial value problems (IVPs) and boundary value problems (BVPs).

We begin our discussion of numerical methods for initial value ordinary differential equations (ODEs) with an introduction of the most basic concepts involved. To illustrate these concepts, we use three simple discretization methods: forward Euler, backward Euler (also called implicit Euler), and trapezoidal. The problem to be solved is written, as before, in the general form

y^′ = f (t, y),    0 ≤ t ≤ b,    (3.1)

with y(0) = c given. You can think of this at first as a scalar ODE—most of what we are going to discuss generalizes to systems directly, and we will highlight occasions where the size of the system is important.

The basic methods developed in Chapter 3 can be adequate for computing approximate solutions of a relatively low accuracy (as we will see, for instance, in Example 4.1) or if the problem being solved is rough in a certain sense (see Section 3.7). But often in practice a quality solution of high accuracy to a relatively smooth problem is sought, and then using a basic, low-order method necessitates taking very small steps in the discretization. This makes the integration process inefficient. Substantially fewer steps are needed when using a higher-order method in such circumstances.

In this chapter we consider another group of methods extending the basic methods of Chapter 3 to higher order; see Figure 4.1 The methods considered here use information from previous integration steps to construct higher-order approximations in a simple fashion. Compared to the Runge-Kutta methods of the previous chapter the methods here typically require fewer function evaluations per step, and they allow a simpler, more streamlined method design, at least from the point of view of order and error estimate. On the other hand, the associated overhead is higher as well, e.g., when wanting to change the step size, and some of the flexibility of one-step methods is lost.

In this chapter and the next two we will consider an ordinary differential equation (ODE) system with m components,

y^′ = f (t, y),   0 < t < b,    (6.1)

subject to m two-point boundary conditions

g(y(0), y(b)) = 0.    (6.2)

We denote the Jacobian matrices of g(u, v) with respect to its first and second argument vectors by

Shooting is a straightforward extension of the initial value techniques that we have seen so far in this book for solving boundary value problems (BVPs). Essentially, one “shoots” trajectories of the same ordinary differential equation (ODE) with different initial values until one “hits” the correct given boundary values at the other interval end. The advantages are conceptual simplicity and the ability to make use of the excellent, widely available, adaptive initial value ODE software. But there are fundamental disadvantages as well, mainly in that the algorithm inherits its stability properties from the stability of the initial value problems (IVPs) that it solves, not just the stability of the given BVP.

As in the previous chapter, we seek numerical methods for boundary value problems (BVPs) based on our knowledge of methods for initial value problems (IVPs). But unlike the previous chapter, here we will not integrate IVPs. Rather, we consider the suitability of the discretizations studied in Chapters 3, 4, and 5 for BVPs. Consider a system of ordinary differential equations (ODEs) of order m,

y^′ = f (t, y),    0 < t < b,    (8.1)

subject to m two-point boundary conditions

g(y(0), y(b)) = 0.    (8.2)

Define a mesh (or a sequence of steps; we refer to the entire mesh as π)

π = {0 = t₀ < t₁ < … < t_N-1 < t_N = b}

with h_n = t_n — t_n-1 the nth step size, and consider solving for

y₀, y₁, …, y_N-1, y_N

with y_n the intended approximation of y(t_n). The following observations are straightforward.

In this chapter and the next we study differential-algebraic equations (DAEs), already introduced in Section 1.3. Here we consider the mathematical structure of such systems and some essential analytical transformations. Numerical approaches and discretizations are discussed in the next chapter. But here, too, our motivation remains finding practical computer solutions. Compared to Chapters 2 and 6, this chapter is unusually long. One reason is that DAE theory is much more recent than ordinary differential equation (ODE) theory. As a result, DAE theory is more in a state of flux, and good expositions are scarce. More importantly, understanding the principles highlighted here is essential for, and will get you a long way towards, constructing good numerical algorithms.

Numerical approaches for the solution of differential-algebraic equations (DAEs) can be divided roughly into two classes: (i) direct discretizations of the given system and (ii) methods which involve a reformulation (e.g., index reduction), combined with a discretization.

The back matter includes the Index