Higher Order Differential Equations
In the preceding two Steps, we have discussed numerical methods for solving the first-order initial value problem y' = f(x, y), y(x0) = y(x0). However, ordinary differential equations which arise in practice are often of higher order. For example, as explained in STEP 1, a more realistic differential equation for the motion of a pendulum is given by
which may have to be solved subject to given values of y(x0) = y0 and y'(x0) = y'0. (In order to ensure notational consistency with the preceding two Steps, we have changed the variables from q and t to y and x, respectively.) In general, an initial value problem for an n-th order differential equation may assume the form:
.
We shall see how this n-th order initial value problem may be rewritten as a system of first-order initial value problems, which leads us to numerical procedures for the solution of the general initial value problem which turn out to be extensions of the numerical methods considered in the preceding two Steps.
If we set wj = y(j-1) for j = 1, 2, . . ., then the n-th order differential equation
becomes y(n)=g(x,y,y',y",···,y(n-2),y(n-1)). Moreover, since w'j = wj+1 for j = 1, 2, . . . , we have an equivalent system of n first-order differential equations, which is subject to the n initial conditions
.
For computational purposes, we choose to regard this hierarchy as a system of n first-order initial value problems. Thus, for example, if the initial conditions for the simple pendulum are y(0) = 0 and y'(0) = 1, then one abtains the system of two first-order initial value problems
.
Note that an initial value problem for a more general system of n first-order differential equations is given by:
for j =1, 2, . . ., n.
For the sake of simplicity, we shall consider only the case n = 2, when the second-order initial value problem
.
yields the first-order system:
.
In order to extend any of the numerical methods, considered in the preceding two Steps, we simply apply the methods to each of these two initial value problems. The easiest way to see how this may be done is to write the system in vector form. Setting:
,
we see that the problem becomes:
.
Now we recall that Euler's method for solving y' = f(x, y) was
.
with the given starting (initial) value y(x0) = y0. The analogous method for solving the system is given by
with w(x0) = w0 given; again, the convention of the preceding steps by which subscripts denote the components xn = x0 + nh has been adopted; moreover, we shall denote the components of wn by wl,n and w2,n. These are the approximations to w1(xn) = y(xn) and w2(xn) = y(xn), respectively. If we write out the separate components in these vectors, we find
.
and
.
As another illustration, we recall the Runge-Kutta method:
.
The analogous method of solving the second-order initial value problem becomes
.
In component form, the equations are:
.
If we use Euler's method to solve the pendulum problem:
.
the resulting equations are
and
.
With w = 1 and h = 0.2, we obtain the values in the table:
.
If we use the Runge-Kutta method, given in the preceding section, we obtain the following values:
Since this Runge-Kutta method is second-order and the Euler method is only first-order, we might expect the values in this table to be more accurate than those displayed in the preceding table. By obtaining very accurate approximations using much smaller values of h, it may be verified that this is indeed the case.
Apply Euler's method with step length h = 0.2 to obtain approximations y(1) and y'(1) for the second-order initial value problem:
.