THE METHOD OF FINITE-DIFFERENCES
HENRY JONES
1. Introduction
In numerical computations of complex differential equations, Finite-Difference
methods in general are used extensively to create sets of mesh points at which
derivatives are approximated and the particular initial or boundary value data are
imposed. Then systems of equations are solved to satisfy the differential equation.
This paper presents a derivation and overview of difference methods for linear and
non-linear, second order ordinary differential equations. With numerical solutions
to partial differential equations being the motivation for this project, discussion
of application of finite-difference methods to the wave equation will serve as a
conclusion.
2. 11.3 Finite-Difference Methods for Linear Problems
In comparison with earlier methods of differential equation approximation in
Chapter 11 of our text, finite-difference methods for linear problems have better sta-
bility characteristics, but require more computation to obtain a specified accuracy.
For boundary-value problems(BVP), finite difference methods replace each deriv-
ative in the ODE with the appropriate difference quotient approximation, many
of which were discussed in chapter 4 including the three and five point- midpoint
formulas for first and second derivatives.
Depending on the equation, the computational equipment, or the physical phe-
nomena being modelled, the difference quotient and step size are chosen to maintain
a certain order of truncation error. Due to the the instability of derivative approxi-
mations resulting from the dominance of round-off error for small step size h values,
our choice of h should not be excessively small.
Discrete Approximation
The finite difference method for the linear, second order BVP,
y
′′
= p(x)y
′
+ q(x)y + r(x) for a ≤ x ≤ b,
with y(a) = α and y(b) = β, necessitates the use of difference-quotient approxima-
tions for y
′
and y
′′
. First we select N such that h = (b − a)/(N + 1) and thus our
mesh points are defined by x
i
= a + ih for i = 0, 1, . . . , N + 1. Therefore at the
mesh points, the differential equation to be approximated is
y
′′
(x
i
) = p(x
i
)y
′
(x
i
) + q(x
i
)y(x
i
) + r(x
i
).
We can then expand y with two, third degree Taylor Polynomials centered at
x
i
and evaluated at x
i+1
and x
i−1
respectively. Note that we must assume that
1
2 HENRY JONES
y ∈ C
4
[x
i−1
, x
i+1
]. For x
i+1
we have
y(x
i+1
) = y(x
i
+ h) = y(x
i
) + hy
′
(x
i
) +
h
2
2
y
′′
(x
i
) +
h
3
6
y
′′′
(x
i
) +
h
4
24
y
(4)
(ξ
+
i
),
for some ξ
+
i
∈ (x
i
, x
i+1
). Similarly we have
y(x
i−1
) = y(x
i
− h) = y(x
i
) − hy
′
(x
i
) +
h
2
2
y
′′
(x
i
) −
h
3
6
y
′′′
(x
i
) +
h
4
24
y
(4)
(ξ
−
i
),
for some ξ
−
i
∈ (x
i−1
, x
i
). If we then add these equations, we see that
y(x
i+1
) + y(x
i−1
) = 2y(x
i
) + h
2
y
′′
(x
i
) +
h
4
24
(y
(4)
(ξ
+
i
) + y
(4)
(ξ
−
i
)),
and thus that
y
′′
(x
i
) =
1
h
2
y(x
i+1
) − 2y(x
i
) + y(x
i−1
)
−
h
4
24
(y
(4)
(ξ
+
i
) + y
(4)
(ξ
−
i
)).
Since we have assumed y ∈ C
4
[x
i−1
, x
i+1
], we can make use of the intermediate
value theorem to say that because y
(4)
∈ C[ξ
−
i
, ξ
+
i
], there is some intermediate
value K between y
(4)
(ξ
−
i
) and y
(4)
(ξ
+
i
) (supposing these are distinct) such that
K =
y
(4)
(ξ
−
i
) + y
(4)
(ξ
+
i
)
2
.
This can be done because the average between two points is always an intermediate
value of the two points. Thus there is some ξ
i
∈ (ξ
−
i
, ξ
+
i
) such that y
(4)
(ξ
i
) = K,
which allows us to simplify our error term to give the centered-difference formula
for y
′′
(x
i
):
y
′′
(x
i
) =
1
h
2
y(x
i+1
) − 2y(x
i
) + y(x
i−1
)
−
h
2
12
y
(4)
(ξ
i
).
From our work in section 4.1, we can obtain a centered-difference formula for
y
′
(x
i
) through a similar derivation with a second order Taylor Polynomial. This
results in
y
′
(x
i
) =
1
2h
[y(x
i+1
) − 2y(x
i
) + y(x
i−1
)] −
h
2
6
y
′′′
(η
i
)
for some η
i
in (x
i−1
, x
i+1
).
Using these approximations for y
′
and y
′′
, our differential equation becomes
y(x
i+1
) − 2y(x
i
) + y(x
i−1
)
h
2
= p(x
i
)
"
y(x
i+1
) − y(x
i−1
)
2h
#
+ q(x
i
)y(x
i
) + r(x
i
)
−
h
2
12
[2p(x
i
)y
′′′
(η
i
) − y
(4)
(ξ
i
)]
at the mesh points for defined by our choice of N . Not that here we have used two,
third order polynomials and thus our error term is of order O(h
2
). Using a fifth
order Taylor Polynomial for this derivation would have given an error term of order
O(h
4
) but is more computationally expensive.
THE METHOD OF FINITE-DIFFERENCES 3
Thus for our approximation scheme, we can define the linear system of equations
given by w
0
= α and w
N+1
= β and
−w
i+1
+ 2w
i
− y
i−1
h
2
+ p(x
i
)
w
i+1
− w
i−1
2h
+ q(x
i
)w
i
= −r(x
i
)
for each i = 1, 2, . . . , N. This equation can be rewritten in the more useful form of
−
1 +
h
2
p(x
i
)
w
i−1
+ (2 + h
2
q(x
i
))w
i
−
1 −
h
2
p(x
i
)
w
i+1
= −h
2
r(x
i
).
Using this form of our ODE approximation, we then define a resulting system of
equations in tridiagonal N × N matrix form Aw = b where
A =
2 + h
2
q(x
1
) −1
h
2
p(x
1
) 0 0
−1
h
2
p(x
2
) 2 + h
2
q(x
2
) −1 +
h
2
p(x
2
)
0 0
−1 +
h
2
p(x
N−1
)
0 0 −1 −
h
2
p(x
N
) 2 + h
2
q(x
N
)
,
w =
w
1
w
2
w
N−1
w
N
, and b =
−h
2
r(x
1
) +
1 +
h
2
p(x
1
)
w
0
−h
2
r(x
2
)
−h
2
r(x
N−1
)
−h
2
r(x
N
) +
1 −
h
2
p(x
N
)
w
N+1
.
This system of linear equations has a unique solution under the conditions of the
following theorem (Theorem 11.3 in the text).
Theorem 1. Suppose the p, q, and r are continuous on [a, b]. If q(x) ≥ 0 on [a, b],
then the tridiagonal linear system above has a unique solution given that h < 2/L
where L = max
a≤x≤b
|p(x)|.
3. Non-Linear BVPs
Having now considered the linear case, we now turn to the case of the general
nonlinear ODE boundary-value problem
y
′′
= f(x, y, y
′
), for a ≤ x ≤ b, with y(a) = α and y(b) = β.
To solve this numeric problem, we must use an iterative process to solve because
the system of equations is no longer guaranteed to be linear.
For the development of the method, we must assume of f, f
x
, and f
′
y
are contin-
uous on the domain where x is bounded and y and y
′
are finite. Additionally, for
some δ > 0, we must have f
y
≥ δ on the domain D and that the max of f
y
and
4 HENRY JONES
f
y
′
are constants on the domain D. By theorem 11.1 in the text, we then know a
unique solution exists.
At the mesh points x
i
, the differential equation is given by
1
h
2
y(x
i+1
)−2y(x
i
)−y(x
i−1
)
!
= f
x
i
, y(x
i
),
1
2h
[y(x
i+1
)−y(x
i−1
)]−
h
2
6
y
′′′
(η
i
)
!
+
h
2
12
y
(4)
(ξ
i
)
for some ξ
i
and η
i
in the interval (x
i−1
, x
i+1
).
Once again the finite difference method arises from omitting the errors terms
and imposing the boundary conditions
w
0
= α, w
N+1
= β,
and solving the new non-linear system of equations
−
w
i+1
− 2w
i
+ w
i−1
h
2
+ f
x
i
, w
i
,
w
i+1
− w
i−1
2h
= 0,
for each i = 1, . . . , N . Provided that h < 2/L, this nonlinear system has a unique
solution which is obtained through Newton’s Method for Iterations detailed in sec-
tion 10.2 of the text. An investigation into this technique is beyond the scope of
this paper.
4. Finite-Difference for the Wave Equation
As a closing to this paper, the previously discussed finite-difference methods for
ODEs will now be extrapolated to the three dimensional, wave equation partial
differential equation given by
∂
2
u
∂t
2
(x, t) − α
2
∂
2
u
∂x
2
(x, t) = 0,
where α is some constant. Following the text’s derivation of the finite-difference
scheme, we impose conditions
0 < x < L, t > 0
u(0, t) = u(l, t) = 0, for t > 0
u(x, 0) = f(x), and
∂u
∂t
(x, 0) = g(x), for 0 ≤ x ≤ l.
Then we select step sizes h and k to create our two-dimensional grid of mesh
points
x
i
= ih and t
j
= jk,
for each i = 0, 1, . . . , m and j = 0, 1, . . .. Thus at the wave equation become
∂
2
u
∂t
2
(x
i
, t
j
) − α
2
∂
2
u
∂x
2
(x
i
, t
j
) = 0.
Following a similar derivation to the linear second order derivative at x
i
, our differ-
ence equation becomes(neglecting error term and summing the centered-difference
quotient for the two partial derivatives)
w
i,j+1
− 2w
i.j
+ w
i,j−1
k
2
− α
2
w
i+1,j
− 2w
i,j
+ w
i−1,j
h
2
= 0.
The details of this method are beyond the scope of this paper, but it is interesting
to note the deep parallels with the centered-difference formulas from our numerical
differentiation methods and connections to the simpler linear ODE derivation of
THE METHOD OF FINITE-DIFFERENCES 5
finite-differences. Ultimately another tridiagonal matrix is a part of our system of
equations which involves more complex matrix matrix multiplication.
5. Conclusion
Despite having initially very high hopes for the end result of this project, I am
still pleased with what I have come to learn and my Python implementation of
algorithm 12.4. I am very interested in physics, the natural sciences, and PDE
theory, so this final project was of great interest to me although I didn’t have
time to take on some of my larger objectives like the shallow water wave equation.
Hopefully I will continue to build my understanding to take on more advanced
projects in this realm in the future.
6. Bibliography
Burden RL, Faires DJ, Burden AM, editors. 2014. Numerical Analysis [Inter-
net]. Boston, MA Cengage Learning;[cited 2022 May 18].
