1 Scope of the Chapter

This chapter is concerned with the numerical solution of ordinary differential equations. There are two main types of problem: those in which all boundary conditions are specified at one point (initial value problems), and those in which the boundary conditions are distributed between two or more points (boundary value problems and eigenvalue problems). Functions are available for initial value problems, two-point boundary value problems and Sturm–Liouville eigenvalue problems.

2 Background to the Problems

For most of the functions in this chapter a system of ordinary differential equations must be written in the form

y1′=f1x,y1,y2,…,yn,y2′=f2x,y1,y2,…,yn,⋮yn′=fnx,y1,y2,…,yn,

that is the system must be given in first-order form. The n dependent variables (also, the solution) y1,y2,…,yn are functions of the independent variable x, and the differential equations give expressions for the first derivatives yi′=dyidx in terms of x and y1,y2,…,yn. For a system of n first-order equations, n associated boundary conditions are usually required to define the solution.

A more general system may contain derivatives of higher order, but such systems can almost always be reduced to the first-order form by introducing new variables. For example, suppose we have the third-order equation

z′′′+zz′′+kl-z′2=0.

We write y1=z, y2=z′, y3=z′′, and the third-order equation may then be written as the system of first-order equations

y1′=y2y2′=y3y3′=-y1y3-kl-y22.

For this system n=3 and we require 3 boundary conditions in order to define the solution. These conditions must specify values of the dependent variables at certain points. For example, we have an initial value problem if the conditions are

y1=0.1 at ​x=0y2=0.1 at ​x=0y3=0.1 at ​x=0.

These conditions would enable us to integrate the equations numerically from the point x=0 to some specified end point. We have a boundary value problem if the conditions are

y1=0 at ​x=0y2=0 at ​x=0y2=1 at ​x=10.

These conditions would be sufficient to define a solution in the range 0≤x≤10, but the problem could not be solved by direct integration (see Section 2.2). More general boundary conditions are permitted in the boundary value case.

2.1 Initial Value Problems

To solve first-order systems, initial values of the dependent variables
yi, for i=1,2,…,n, must be supplied at a given point, a. Also a point, b, at which the values of the dependent variables are required, must be specified. The numerical solution is then obtained by a step-by-step calculation which approximates values of the variables
yi, for i=1,2,…,n, at finite intervals over the required range a,b. The functions in this chapter adjust the step length automatically to meet specified accuracy tolerances. Although the accuracy tests used are reliable over each step individually, in general an accuracy requirement cannot be guaranteed over a long range. For many problems there may be no serious accumulation of error, but for unstable systems small perturbations of the solution will often lead to rapid divergence of the calculated values from the true values. A simple check for stability is to carry out trial calculations with different tolerances; if the results differ appreciably the system is probably unstable. Over a short range, the difficulty may possibly be overcome by taking sufficiently small tolerances, but over a long range it may be better to try to reformulate the problem.

A special class of initial value problems are those for which the solutions contain rapidly decaying transient terms. Such problems are called stiff; an alternative way of describing them is to say that certain eigenvalues of the Jacobian matrix ∂fi∂yj have large negative real parts when compared to others. These problems require special methods for efficient numerical solution; the methods designed for non-stiff problems when applied to stiff problems tend to be very slow, because they need small step lengths to avoid numerical instability. A full discussion is given in Hall and Watt (1976) and a discussion of the methods for stiff problems is given in Berzins et al. (1988).

2.2.1 Finite difference methods

If a boundary value problem seems insoluble by the above methods and a good estimate for the solution of the problem is known at all points of the range then a finite difference method may be used. Finite difference equations are set up on a mesh of points and estimated values for the solution at the grid points are chosen. Using these estimated values as starting values a Newton iteration is used to solve the finite difference equations. The accuracy of the solution is then improved by deferred corrections or the addition of points to the mesh or a combination of both. The method does not suffer from the difficulties associated with the shooting method but good initial estimates of the solution may be required in some cases and the method is unlikely to be successful when the solution varies very rapidly over short ranges. A discussion is given in Chapters 9 and 11 of Gladwell and Sayers (1980) and Chapter 4 of Gladwell (1979a).

3 Recommendations on Choice and Use of Available Functions

There are no functions which deal directly with complex equations. These may however be transformed to larger systems of real equations of the required form. Split each equation into its real and imaginary parts and solve for the real and imaginary parts of each component of the solution. Whilst this process doubles the size of the system and may not always be appropriate it does make available for use the full range of functions provided presently.

There is a simple driving function nag_ode_ivp_adams_gen (d02cjc), which integrates a system over a range and, optionally, computes intermediate output and/or determines the position where a specified function of the solution is zero.

3.1.3 BDF functions

General functions for implicit ordinary differential equations with options for forms of numerical linear algebra are provided by a suite of functions using the DASSL implementation (see Brenan et al. (1996)) of the BDF. The main solver in this suite is nag_dae_ivp_dassl_gen (d02nec) is designed for solving systems of the form,

Ft,y,y′=0.

These formulations permits solution of differential/algebraic systems (DAEs). Additionally nag_dae_ivp_dassl_gen (d02nec) can be used to solve difficult algebraic problems by continuation; for example, the nonlinear algebraic problem

General functions for explicit and implicit ordinary differential equations with a wide range of options for integrator choice and special forms of numerical linear algebra are provided in sub-chapter d02m–n. A separate document describing the use of this sub-chapter is given immediately before the functions of the sub-chapter.

There is a simple driving function nag_ode_ivp_bdf_gen (d02ejc), which integrates a system over a range and, optionally, computes intermediate output and/or determines the position where a specified function of the solution is zero. It has a specification similar to the Adams function nag_ode_ivp_adams_gen (d02cjc) except that to solve the equations arising in the BDF method an approximation to the Jacobian ∂fi∂yj is required. This approximation can be calculated internally but you may supply an analytic expression. In most cases supplying a correct analytic expression will reduce the amount of computer time used.

3.2 Boundary Value Problems

For simple boundary value problems with assigned boundary values you may prefer to use a code based on the finite difference method for which there is a function with simple calling sequence (nag_ode_bvp_fd_nonlin_fixedbc (d02gac)).

For difficult boundary value problems, where you need to exercise some control over the calculation, and where the collocation method proves unsuccessful, you may wish to try the alternative
method of
finite differences (nag_ode_bvp_fd_nonlin_gen (d02rac)).

Note that it is not possible to make a fully automatic boundary value function, and you should be prepared to experiment with different starting values or a different function if the problem is at all difficult.

3.2.1 Finite difference methods

You may find that convergence is difficult to achieve using nag_ode_bvp_fd_nonlin_fixedbc (d02gac) since only specifying the unknown boundary values and the position of the finite difference mesh is permitted. In such cases you may use nag_ode_bvp_fd_nonlin_gen (d02rac), which permits specification of an initial estimate for the solution at all mesh points and allows the calculation to be influenced in other ways too. nag_ode_bvp_fd_nonlin_gen (d02rac) is designed to solve a general nonlinear two-point boundary value problem with nonlinear boundary conditions.

A function, nag_ode_bvp_fd_lin_gen (d02gbc), is also supplied specifically for the general linear two-point boundary value problem written in a standard ‘textbook’ form.

You are advised to use interpolation functions from Chapter e01 to obtain solution values at points not on the final mesh.

3.2.2 Chebyshev integration method

The Chebyshev integration method is an implementation of the Chebyshev collocation method which is fully described and compared against other implementations in Muite (2010). nag_ode_bvp_ps_lin_solve (d02uec) solves a linear constant coefficient boundary value problem using the Chebyshev integration formulation on a Chebyshev Gauss–Lobatto grid and solving in the coefficient space. The required Chebyshev Gauss–Lobatto grid points on a given arbitrary interval a,b can first be generated using nag_ode_bvp_ps_lin_cgl_grid (d02ucc). Then nag_ode_bvp_ps_lin_coeffs (d02uac) obtains the Chebyshev coefficients for the right-hand side (of system) function discretized on the obtained Chebyshev Gauss–Lobatto grid. nag_ode_bvp_ps_lin_solve (d02uec) then solves the problem in Chebyshev coefficient space using the integration formulation. Finally nag_ode_bvp_ps_lin_cgl_vals (d02ubc) evaluates the solution (or one of its lower order derivatives) from the set of Chebyshev coefficients returned by nag_ode_bvp_ps_lin_solve (d02uec) on the Chebyshev Gauss–Lobatto grid on a,b. The set of functions can be used to solve up to fourth order boundary value problems.

3.3 Summary of Recommended Functions

Problem

Function

RK Method

Adams Method

BDF Method

Initial Value ProblemsDriver Functions

Integration over a range with optional intermediate output and optional determination of position where a function of the solution becomes zero