13/04/2026, 16:50 Finite Difference Method — Python Numerical Methods
This notebook contains an excerpt from the Python Programming and Numerical Methods - A Guide
for Engineers and Scientists, the content is also available at Berkeley Python Numerical Methods.
The copyright of the book belongs to Elsevier. We also have this interactive book online for a better
learning experience. The code is released under the MIT license. If you find this content useful, please
consider supporting the work on Elsevier or Amazon!
< 23.2 The Shooting Method | Contents | 23.4 Numerical Error and Instability >
Finite Difference Method
Another way to solve the ODE boundary value problems is the finite difference method, where we can use finite
difference formulas at evenly spaced grid points to approximate the differential equations. This way, we can
transform a differential equation into a system of algebraic equations to solve.
In the finite difference method, the derivatives in the differential equation are approximated using the finite
difference formulas. We can divide the the interval of [a, b] into n equal subintervals of length h as shown in the
following figure.
Commonly, we usually use the central difference formulas in the finite difference methods due to the fact that
they yield better accuracy. The differential equation is enforced only at the grid points, and the first and second
derivatives are:
dy y i+1 − y i−1
=
dx 2h
2
d y y i−1 − 2y i + y i+1
=
2 2
dx h
These finite difference expressions are used to replace the derivatives of y in the differential equation which leads
to a system of n + 1 linear algebraic equations if the differential equation is linear. If the differential equation is
nonlinear, the algebraic equations will also be nonlinear.
EXAMPLE: Solve the rocket problem in the previous section using the finite difference method, plot the altitude of
the rocket after launching. The ODE is
2
d y
= −g
dt2
with the boundary conditions y(0) = 0 and y(5) = 50 . Let’s take n = 10 .
Since the time interval is [0, 5] and we have n = 10 , therefore, h = 0.5 , using the finite difference approximated
derivatives, we have
y0 = 0
[Link] the finite difference method%2C,using the finite difference formulas. 1/5
13/04/2026, 16:50 Finite Difference Method — Python Numerical Methods
2
y i−1 − 2y i + y i+1 = −gh , i = 1, 2, . . . , n − 1
y 10 = 50
if we use matrix notation, we will have:
1 0
⎡ ⎤ y0 0
⎡ ⎤ ⎡ ⎤
⎢ 1 −2 1 ⎥ 2
⎢ ⎥ ⎢ y1 ⎥ ⎢ −gh ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢ ⎥ = ⎢ ⎥
⎢
⋱ ⋱ ⋱
⎥ ⎢ ... ⎥ ⎢ ... ⎥
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥
⎢y ⎥ ⎢ −gh2 ⎥
⎢ 1 −2 1⎥ n−1
⎣ ⎦ ⎣ ⎦
⎣ ⎦ yn 50
1
Therefore, we have 11 equations in the system, we can solve it using the method we learned in chapter 14.
import numpy as np
import [Link] as plt
[Link]('seaborn-poster')
%matplotlib inline
n = 10
h = (5-0) / n
# Get A
A = [Link]((n+1, n+1))
A[0, 0] = 1
A[n, n] = 1
for i in range(1, n): Print to PDF
A[i, i-1] = 1
A[i, i] = -2
A[i, i+1] = 1
print(A)
# Get b
b = [Link](n+1)
b[1:-1] = -9.8*h**2
b[-1] = 50
print(b)
# solve the linear equations
y = [Link](A, b)
t = [Link](0, 5, 11)
[Link](figsize=(10,8))
[Link](t, y)
[Link](5, 50, 'ro')
[Link]('time (s)')
[Link]('altitude (m)')
[Link]()
[[ 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 1. -2. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 1. -2. 1. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 1. -2. 1. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 1. -2. 1. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 1. -2. 1. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 1. -2. 1. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 1. -2. 1. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 1. -2. 1. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 1. -2. 1.]
[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]]
[ 0. -2.45 -2.45 -2.45 -2.45 -2.45 -2.45 -2.45 -2.45 -2.45 50. ]
[Link] the finite difference method%2C,using the finite difference formulas. 2/5
13/04/2026, 16:50 Finite Difference Method — Python Numerical Methods
dy yi+1 −yi−1
Now, let’s solve y ′ (0), from the finite difference formula, we know that dx
=
2h
, which means that
y1 −y−1
′
y (0) =
2h
, but we don’t know what is y −1 . Actually, we can calculate y −1 since we know the y values on
y−1 −2y0 +y1
each grid point. From the 2nd derivative finite difference formula, we know that 2
= −g , therefore, we
h
can solve for y −1 and then get the launching velocity. See the calculation below.
y_n1 = -9.8*h**2 + 2*y[0] - y[1]
(y[1] - y_n1) / (2*h)
34.5
We can see that we get the correct launching velocity using the finite difference method. To make you more
comfortable with the method, let’s see another example.
TRY IT! Using finite difference method to solve the following linear boundary value problem
′′
y = −4y + 4x
with the boundary conditions as y(0) = 0 and y ′ (π/2) = 0 . The exact solution of the problem is
y = x − sin2x , plot the errors against the n grid points (n from 3 to 100) for the boundary point y(π/2) .
Using the finite difference approximated derivatives, we have
y0 = 0
2
y i−1 − 2y i + y i+1 − h (−4y i + 4xi ) = 0, i = 1, 2, . . . , n − 1
2
2y n−1 − 2y n − h (−4y n + 4xn ) = 0
yn+1 −yn−1
The last equation is derived from the fact that = 0 (the boundary condition y ′ (π/2) = 0 ). Therefore,
2h
y n+1 = y n−1 .
if we use matrix notation, we will have:
1 0
⎡ ⎤ y0 0
⎡ ⎤ ⎡ ⎤
2
⎢ 1 −2 + 4h 1 ⎥ 2
⎢ ⎥ ⎢ y1 ⎥ ⎢ 4h x1 ⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥ = ⎢ ⎥
⋱ ⋱ ⋱ ⎢ ... ⎥ ⎢
...
⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢y ⎥ ⎢ 4h2 x ⎥
2
⎢ 1 −2 + 4h 1 ⎥ n−1 n−1
⎣ ⎦ ⎣ 2 ⎦
⎣ 2 ⎦ yn 4h xn
2 −2 + 4h
[Link] the finite difference method%2C,using the finite difference formulas. 3/5
13/04/2026, 16:50 Finite Difference Method — Python Numerical Methods
def get_a_b(n):
h = ([Link]/2-0) / n
x = [Link](0, [Link]/2, n+1)
# Get A
A = [Link]((n+1, n+1))
A[0, 0] = 1
A[n, n] = -2+4*h**2
A[n, n-1] = 2
for i in range(1, n):
A[i, i-1] = 1
A[i, i] = -2+4*h**2
A[i, i+1] = 1
# Get b
b = [Link](n+1)
for i in range(1, n+1):
b[i] = 4*h**2*x[i]
return x, A, b
x = [Link]/2
v = x - [Link](2*x)
n_s = []
errors = []
for n in range(3, 100, 5):
x, A, b = get_a_b(n)
y = [Link](A, b)
n_s.append(n)
e = v - y[-1]
[Link](e)
[Link](figsize = (10,8))
[Link](n_s, errors)
[Link]('log')
[Link]('n gird points')
[Link]('errors at x = $\pi/2$')
[Link]()
We can see with denser grid points, we are approaching the exact solution on the boundary point.
The finite difference method can be also applied to higher-order ODEs, but it needs approximation of the higher-
order derivatives using the finite difference formula. For example, if we are solving a fourth-order ODE, we will
need to use the following:
4
d y y i−2 − 4y i−1 + 6y i − 4y i+1 + y i+2
=
4 4
dx h
We won’t talk more on the higher-order ODEs, since the idea behind to solve it is similar to the second-order ODE
we discussed above.
< 23.2 The Shooting Method | Contents | 23.4 Numerical Error and Instability >
[Link] the finite difference method%2C,using the finite difference formulas. 4/5
13/04/2026, 16:50 Finite Difference Method — Python Numerical Methods
© Copyright 2020.
[Link] the finite difference method%2C,using the finite difference formulas. 5/5