CHAPTER

8
NUMERICAL SOLUTION OF ORDINARY
DIFFERENTIAL EQUATIONS

8.1 INTRODUCTION
Numerical methods are becoming more and more important in engineering
applications, simply because of the difficulties encountered in finding exact
analytical solutions but also, because of the ease with which numerical tech-
niques can be used in conjunction with modern high-speed digital computers.
Several numerical procedures for solving initial value problems involving first-
order ordinary differential equations are discussed in this chapter.
In spite of the fact that the error analysis is an important part of any
numerical procedure, the discussion in this chapter is limited primarily to the
use of the procedure itself. The theory of errors and error analysis is some-
times fairly complex and goes beyond the intended scope of this chapter.
An ordinary differential equation is one in which an ordinary derivative of
a dependent variable y with respect to an independent variable x is related in
a prescribed manner to x, y, and lower derivatives. The most general form of
an ordinary differential equation of nth order is given by

dn y  dy d 2 y d n −1 y 
= f  x, y , , 2 ,...., n −1  (8.1)
dx n
 dx dx dx 

The Equation (8.1) is termed ordinary because there is only one inde-
pendent variable.

ANMF.CH08_4PP.indd 557 3/1/2023 2:36:14 PM

 558 • APPLIED NUMERICAL METHODS USING MATLAB

To solve an equation of the type (Equation (8.1)), we also require a set

of conditions. When all the conditions are given at one value x and the solu-
tion proceeds from that value of x, we have an initial-value problem. When
the conditions are given at different values of x, we have a boundary-value
problem.
A general solution of an ordinary differential equation (Equation (8.1))
would be a relation between y, x, and n arbitrary constants which is of form

f(x, y, c1, c2, …, cn) = 0 (8.2)

If particular values are given to the constants cn in Equation (8.2), then

the resulting solution is called a particular solution. There are many analyti-
cal methods available for finding the solution of the Equation (8.1). However,
there exist a large number of ordinary differential equations in science and
engineering, whose solutions cannot easily be obtained by the well-known
analytical methods. For such ordinary differential equations, we can obtain an
approximate solution of a given ordinary differential equations using numeri-
cal methods under the given initial conditions.
Any ordinary differential equation can be replaced by a system of first-
order differential equations (which involve only first derivatives). The single
first-order ordinary differential equation with an initial value is a special case
of Equation (8.1). It is described by
dy
= f (x , y ) y = y0 at x = x0 (8.3)
dx
The description in Equation (8.3) consists of the differential equation
itself and a given solution y0 at initial location x0. We then obtain the solution
y as x ranges from its initial value to some other value.
The general solution of Equation (8.3) can be obtained in two forms:

1. the values of y as a power series in independent variable x

2. as a set of tabulated values of x and y.
There are two categories of methods to solve ordinary differential
equations:

1. One-step methods or single-step methods.

2. Step-by-step methods or marching methods.

ANMF.CH08_4PP.indd 558 3/1/2023 2:36:17 PM

 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS • 559

In one-step methods or single-step methods, the information about the

curve represented by an ordinary differential equation at one point is uti-
lized and the solution is not iterated. In step-by-step methods or the march-
ing methods, the next point on the curve is evaluated in short steps ahead, for
equal intervals of width h of the independent variable, by performing itera-
tions till the desired level of accuracy is obtained.
In general, we divide the interval (a, b) on which the solution is derived
into a finite number of subintervals by the points a = x0 < x1 < x2, … < xn = b,
called the mesh points. This is done by setting up xn = x0 + nh.
The existence of the uniqueness of the solution to an initial value problem
in (x0, b) is based on Lipschitz theorem. Lipschitz theorem states that:
a) If f(x, y) is a real function defined and continuous in (x0, b), y∈(−∞, +∞),
where x0 and b are finite.
b) There exists a constant k > 0 called Lipschitz constant such that for any
two values y = y1 and y = y2

|f(x, y1) − (f(x, y2)| < k|k1 − k2|

where x∈(x0, b), then for any y(x0) = y0, the initial value problem [Equation (8.3)],
has unique solution for x∈(x0, b).
Also, there are two types of methods, explicit and implicit, can be used to
compute the solution at each step. Explicit methods are those methods that
use an explicit formula for calculating the value of the dependent variable at
the next value of the independent variable. In an explicit method, the right-
hand side of the equation only has all known quantities. Therefore, the next
unknown value of the dependent variable, yn+1, is calculated by evaluating an
expression of the form:

yn+1 = F(xn, xn+1, yn) (8.4)

where xn, yn and xn+1 are all known quantities.

In implicit methods, the equation used for computing yn+1 from the known
xn, yn and yn+1 has the form:

yn+1 = F(xn, xn+1, yn+1) (8.5)

Here, the unknown yn+1 appears on both sides of the equation. Generally
speaking, the right-hand side of Equation (8.3c) is nonlinear. Therefore, the

ANMF.CH08_4PP.indd 559 3/1/2023 2:36:18 PM

 560 • APPLIED NUMERICAL METHODS USING MATLAB

Equation (8.5) must be solved for yn+1 using suitable numerical methods. In
general, implicit methods give better accuracy over explicit methods at the
expense of additional effort.
In this chapter, one-step or single-step methods such as Picard’s method
of successive approximations and the Taylor series methods were presented.
Euler’s method, the modified Euler’s method, and the Runge-Kutta methods
of order two and four were also presented. The Adam-Moulton predictor-
corrector method and Milne’s predictor-corrector methods were presented
among the step-by-step methods or the marching methods. All of these meth-
ods will be illustrated with worked examples.

8.2 ONE-STEP METHODS OR SINGLE-STEP METHODS

In single-step explicit method, the approximate solution (xn+1, yn+1) is com-
puted from the known solution at point (xn, yn) using
xn+1 = xn + h (8.6)
yn+1 = yn + (slope) h (8.7)
This is illustrated in Figure 8.1. Here in Equation (8.6), h is the step size
dy
and the slope is a constant that estimates the value of in the interval from
dx
xn to xn+1. The numerical solution starts at the point where the initial value is
known corresponding to n = 1 and point (x1, y1). Then, n is increased to n = 2,
and the solution at the next point, (x2, y2) is computed using Equations (8.6)
and (8.7). This procedure is repeated for n = 3 and so on until the points cover
the whole domain of the solution.

y
y(x)

Exact solution

yn+1
Slope
Numerical
solution

yn

x
xn h xn+1

FIGURE 8.1 Single-step explicit methods.

ANMF.CH08_4PP.indd 560 3/1/2023 2:36:19 PM