Scientific Machine Learning:
Numerical Analysis
Physical laws are written as derivatives, integrals and differential equations, and a computer can handle none of these directly. Numerical analysis supplies the approximations that make them computable: roots of nonlinear equations, finite difference derivatives, numerical integration, and time-stepping schemes for ordinary and partial differential equations. These are the classical tools that scientific machine learning both competes with and builds on, and every one of them reappears later, inside a loss function or as the reference solution a network is compared against.
Learning Objectives¶
After studying this notebook, you should be able to
Roots of nonlinear equations
- state the rootfinding problem $f(r) = 0$, and explain why the residual $f(\tilde r)$, not the error $\tilde r - r$, is the quantity a numerical method can actually measure;
- apply the bisection method, write down its error bound after $k$ steps, and explain why it always converges and why it is nevertheless slow;
- run a fixed-point iteration $x_{k+1} = g(x_k)$ and predict from $|g'(p)|$ whether it converges, and how fast;
- derive Newton's method from the tangent line, implement it with sensible stopping criteria, and recognize its quadratic convergence in a table of errors;
- derive the secant method, explain why it needs no derivative, and compare its cost per step with Newton's;
- compare the three methods on the same equation in terms of iterations and function evaluations, and choose between them.
Finite difference derivatives
- explain why a derivative must be approximated on a computer, and why the limit $\Delta x \to 0$ cannot simply be taken;
- derive the forward, backward and central difference formulas from Taylor series, and read off the order of accuracy of each from the leading error term;
- derive the central formula for the second derivative, and build a higher-order stencil by combining two lower-order ones;
- measure an order of accuracy numerically from the slope of an error curve on a log-log plot;
- differentiate a function sampled on a grid, and handle the two boundaries with one-sided formulas.
Numerical integration
- explain why an integral must often be computed numerically, and recognize the Riemann sum hidden in the definition of the integral;
- write down the left-sided and right-sided rectangle rules, implement both, and explain why they are first-order accurate;
- obtain the trapezoidal rule as the average of the two rectangle rules, and explain why averaging raises the order from one to two;
- estimate an area by hit-or-miss sampling and an integral by the sample mean, and explain why a uniform random sample turns an integral into an expected value;
- state the $O(N^{-1/2})$ convergence rate of a Monte Carlo estimate, explain why it does not depend on the dimension, and say at which dimension it overtakes a grid-based rule.
Numerical solutions to ODEs
- write a differential equation as a first-order system $\dot{x} = f(x)$ with an initial condition, and explain what a numerical solver actually produces;
- derive the forward and backward Euler methods from the forward and backward difference formulas, and write both in matrix form for a linear system;
- implement both methods as a loop and compare them with a library solver;
- explain the difference between local and global error, and show why both Euler methods are first-order accurate;
- determine, from $|1 + \lambda\Delta t|$ or from the eigenvalues of $I + \Delta t A$, whether a step size is stable, and explain why the implicit method has no such restriction;
- set up and simulate a nonlinear system such as the pendulum, and recognize the signature of an unstable integrator in the result.
Boundary value problems
- explain why a boundary value problem cannot be marched like an initial value problem, and why its unknowns must be solved for simultaneously;
- discretize $u'' = f$ with the three-point central stencil and assemble the resulting tridiagonal system;
- move the known boundary values to the right-hand side, solve the system, and confirm second-order accuracy by grid refinement;
- state how a Neumann condition changes the system, and why the order of the boundary formula matters.
Partial differential equations
- classify a second-order PDE as hyperbolic, parabolic or elliptic, and say which physical behaviour each describes;
- derive the explicit FTCS scheme for the heat equation by combining a forward difference in time with a central difference in space;
- state the stability condition $r = \alpha\Delta t/\Delta x^2 \le 1/2$, explain it both from the sign of the weight $1 - 2r$ and from the amplification factor, and recognize the sawtooth that appears when it is violated;
- derive the implicit BTCS scheme, assemble the tridiagonal system it requires at every time step, and explain why it is unconditionally stable;
- choose between an explicit and an implicit method for a given problem.
Setup
NumPy, Matplotlib and the SciPy routines used for reference solutions.
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import root_scalar
from scipy.integrate import solve_ivp
from scipy.linalg import expm
1.1.1. Definition¶
Given a continuous function $f$ of one variable, the root finding problem is to find a real number $r$, called a root, such that
$$ f(r) = 0 . $$
Any equation of the form $g(x) = h(x)$ can be written in this standard form by setting $f = g - h$, so the definition is not restrictive.
Roots of a general $f$ cannot be computed exactly in finitely many arithmetic operations. Every method of this chapter is therefore iterative: it produces a sequence $x_1, x_2, x_3, \dots$ that (we hope) converges to $r$, and we stop when the sequence has settled down to the accuracy we need.
1.1.2. A Running Example¶
The methods of Sections 1.2, 1.4 and 1.5 are all compared on one equation,
$$ f(x) = x e^{x} - 2 = 0 , $$
so that they can be judged on equal terms. It cannot be solved in closed form with elementary functions, but a glance at the graph shows that $f$ is increasing on $[0, 1]$ with $f(0) = -2 < 0$ and $f(1) = e - 2 \approx 0.718 > 0$, so there is exactly one root in that interval.
The first code cell imports everything used below and plots $f$. The accurate reference value $r$ printed there is obtained with SciPy's root_scalar; we will use it only to measure the error of our own iterations, which is a luxury unavailable in a real computation.
def f(x):
return x*np.exp(x) - 2
# accurate reference value, used only to measure the error of our own iterations
r = root_scalar(f, bracket = [0, 1], xtol = 1e-15).root
print('reference root r = %.5f' % r)
print('residual f(r) = %.1e' % f(r))
# Plot
xx = np.linspace(0, 1.5, 400)
plt.figure(figsize = (6, 4))
plt.plot(xx, f(xx), 'k', linewidth = 2)
plt.plot([0, 1.5], [0, 0], 'k', linewidth = 0.5)
plt.plot(r, 0, 'ro', markersize = 8)
plt.grid(alpha = 0.2)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('$f(x) = xe^x - 2$')
plt.show()
1.2.1. The Idea¶
The bisection method rests on the intermediate value theorem: if $f$ is continuous on $[a, b]$ and
$$ f(a)\,f(b) < 0 , $$
then $f$ changes sign on the interval, so $[a, b]$ contains at least one root. Such an interval is called a bracket.
The bisection method shrinks the bracket by halving it. Let $x_k$ be the midpoint of the current interval. Evaluating $f(x_k)$ tells us which half still has the sign change: if $f(a)$ and $f(x_k)$ have opposite signs the root lies in $[a, x_k]$, otherwise it lies in $[x_k, b]$. Replace the interval by that half and repeat. The bracket is never lost, so the method cannot fail.
The implementation below returns the sequence of midpoints, so that the convergence can be inspected afterwards. Each pass costs exactly one new evaluation of $f$.
a = 0
b = 1
print('f(a) = %8.4f' % f(a))
print('f(b) = %8.4f' % f(b))
# Store initial bracket for bound calculation
a_initial = a
b_initial = b
# Initial setup
fa = f(a)
xb = [] # List to store history of midpoints
# Halve the bracket until it is narrow enough
while (b - a) / 2 > 1e-14:
xm = (a + b) / 2
fm = f(xm)
xb.append(xm)
# Check which half contains the sign change
if np.sign(fm) == np.sign(fa):
a = xm
fa = fm # Update f(a) for the next iteration
else:
b = xm
xb = np.array(xb)
errb = np.abs(xb - r)
# Display convergence history
print('%3s %15s %12s %12s' % ('k', 'x_k', '|error|', 'bound'))
print('-' * 46)
for k in range(min(10, len(xb))):
bound = (b_initial - a_initial) / (2.0 ** (k + 1)) # 2^-(k+1) for unit interval
print('%3d %15.12f %12.2e %12.2e' % (k + 1, xb[k], errb[k], bound))
# Only print final error if xb is not empty
if len(xb) > 0:
print('\nTotal steps: %s, Final error: %.2e' % (len(xb), errb[-1]))
else:
print('\nNo steps were taken as the initial bracket was already narrow enough.')
The first few midpoints are nowhere near the root — the method knows only the sign of $f$, never how close to zero it is — and the error column is not even monotone: step $4$ is worse than step $3$, and step $8$ worse than step $7$. What is monotone is the last column, the width of the bracket, which halves without fail. Every line of the table is therefore a guarantee rather than a hope, and the next subsection turns that guarantee into a formula.
1.2.2. Convergence¶
The initial bracket has width $b - a$; after $k$ halvings it has width $(b-a)/2^{k}$, and the midpoint $x_k$ is within half of that of the root:
$$ |x_k - r| \le \frac{b - a}{2^{k}} . $$
This is an error bound available before the computation starts. Solving for $k$ gives the number of steps needed for a required accuracy $\delta$,
$$ k \ge \log_2\!\frac{b - a}{\delta} , $$
independently of $f$. For $[a, b] = [0, 1]$ and $\delta = 10^{-14}$ this is $k \approx 46$ — and the loop above indeed stops after $46$ steps, with a final error of about $5 \times 10^{-15}$.
The error is (at worst) multiplied by $\tfrac12$ at every step, which is linear convergence with rate $\sigma = 1/2$: a straight line on a semilogarithmic plot, and about $\log_{10} 2^{-1} \approx 0.3$ decimal digits gained per iteration, i.e. one digit every $3.3$ evaluations of $f$. That rate is fixed: it can be beaten, but only by a method that looks at the values of $f$ and not merely at its sign.
k = np.arange(1, len(xb) + 1)
plt.figure(figsize = (6, 4))
plt.semilogy(k, errb, 'ko-', linewidth = 1.2, markersize = 3)
plt.semilogy(k, 2.0**(-k), 'r--', linewidth = 1.2)
plt.legend(['bisection error', 'bound $2^{-k}$'])
plt.grid(True)
plt.xlabel('$k$')
plt.ylabel('$|x_k - r|$')
plt.show()
Remark (what bisection cannot do). It needs a bracket to start, and finding one requires some knowledge of the function — a plot, physical reasoning, or a scan over a grid. It also cannot find a root at which $f$ touches zero without crossing it, such as $x = 0$ for $f(x) = x^2$, because no bracket exists there. Finally, it uses only the sign of $f$, and therefore never converges faster than rate $1/2$ no matter how smooth or well-behaved the function is. The methods that follow use the values of $f$, gaining speed and losing the guarantee.
1.3.1. From Root Finding to a Fixed-Point Problem¶
We want to solve
$$ f(x) = 0 . $$
Instead of attacking this directly, we rewrite it as a fixed-point problem: find a point that a function maps to itself,
$$ x = g(x) . $$
A point $r$ with $r = g(r)$ is called a fixed point of $g$. The rewriting is always possible algebraically, for instance
$$ \begin{align*} x &= x + f(x) = g(x), \quad \text{or} \\ x &= x - f(x) = g(x) . \end{align*} $$
The choice of $g$ is not unique, and it is not innocent. Every valid rewriting has the same solution, but as the examples below show, some of them produce an iteration that converges quickly, some converge slowly, and some do not converge at all. Choosing $g$ well is the entire craft of this method.
The algorithm
- Choose an initial guess $x_0$.
- Repeat $x_{k+1} = g(x_k)$ until the change between successive iterates is small enough.
That is the whole method: evaluate the function, feed the answer back in, repeat. No derivatives, no linear systems, no bracketing.
1.3.2. When Does It Converge?¶
Let $r$ be the exact solution, so that $r = g(r)$, and define the error at step $k$ as $e_k = x_k - r$. Subtracting the two relations,
$$ \begin{aligned} x_{k+1} &= g(x_k) \\ r &= g(r) \\[4pt] e_{k+1} = x_{k+1} - r &= g(x_k) - g(r) . \end{aligned} $$
By the mean value theorem there is a point $\eta$ between $x_k$ and $r$ with
$$ g(x_k) - g(r) = g'(\eta)\,(x_k - r), $$
so the errors obey the simple recursion
$$ e_{k+1} = g'(\eta)\, e_k \qquad\Longrightarrow\qquad \lvert e_{k+1}\rvert \le \lvert g'(\eta)\rvert \cdot \lvert e_k \rvert . $$
Therefore, if $\lvert g'\rvert < 1$ near the root, every step shrinks the error by roughly that factor and the iteration converges. If $\lvert g'\rvert > 1$, every step enlarges the error and the iteration diverges.
$$ \text{converges if } \lvert g'(r)\rvert < 1, \qquad \text{diverges if } \lvert g'(r)\rvert > 1\ $$
Three further readings of this result are worth keeping in mind.
The convergence is linear: the error is multiplied by approximately the constant $\lvert g'(r)\rvert$ at every step, so the number of correct digits grows by a fixed amount per iteration. A smaller $\lvert g'(r)\rvert$ means faster convergence.
The special case $g'(r) = 0$ is much better than linear. The first-order term vanishes, the next term in the Taylor expansion takes over, and the error is then squared at each step. Section 1.3.4 contains an example.
The result is local. It describes the behaviour once the iterate is near $r$, where $g'(\eta)\approx g'(r)$. A poor starting point can still send the iteration somewhere else entirely.
The geometric picture
The condition has a picture that is worth more than the algebra. Plot $y = g(x)$ and the diagonal $y = x$ on the same axes. A fixed point is where the two curves cross.
Starting at $x_k$ on the horizontal axis, one iteration does two things: move vertically to the curve, which computes $g(x_k)$, and then move horizontally to the diagonal, which turns that value into the next input $x_{k+1}$. Repeating draws a staircase or a spiral, and whether it winds inward toward the crossing or outward away from it is decided by how steep $g$ is there compared with the diagonal. This is the cobweb diagram, and we draw it for every example below.
import numpy as np
import matplotlib.pyplot as plt
A helper for the cobweb diagram. Given a function $g$, a starting point, and a number of steps, it plots $g$, the diagonal, and the alternating vertical and horizontal moves of the iteration.
def cobweb(g, x0, n_steps, a, b, name):
x = np.linspace(a, b, 200)
plt.plot(x, g(x), 'b', linewidth = 2, label = name)
plt.plot(x, x, 'k', linewidth = 1, label = 'y = x')
xk = x0
for k in range(n_steps):
xn = g(xk)
plt.plot([xk, xk], [xk, xn], 'r', linewidth = 1) # vertical: evaluate g
plt.plot([xk, xn], [xn, xn], 'r', linewidth = 1) # horizontal: feed it back
xk = xn
plt.plot(x0, x0, 'ro')
plt.xlabel('x'); plt.ylabel('y'); plt.legend();
1.3.3. Example: Solving $x = \cos(x)$¶
The equation $\cos(x) = x$ cannot be solved in closed form, but it is already in fixed-point form with $g(x) = \cos(x)$. First look at where the solution is.
x = np.linspace(-1.5*np.pi, 1.5*np.pi, 100)
plt.figure(figsize = (6, 4))
plt.plot(x, np.cos(x), linewidth = 2, label = 'cos(x)')
plt.plot(x, x, linewidth = 2, label = 'x')
plt.axvline(x = 0, color = 'k', linestyle = '--')
plt.axhline(y = 0, color = 'k', linestyle = '--')
plt.axis('equal')
plt.ylim([-4, 4]);
plt.xlabel('x');
plt.legend();
plt.show()
The two curves cross exactly once, somewhere near $x = 0.74$. There is one solution and we now compute it.
A first attempt. The most direct way to express "apply $\cos$ again and again" is to write it out by hand.
x = 0.3
print(np.cos(x))
print(np.cos(np.cos(x)))
print(np.cos(np.cos(np.cos(x))))
print(np.cos(np.cos(np.cos(np.cos(x)))))
print(np.cos(np.cos(np.cos(np.cos(np.cos(x))))))
The numbers are clearly settling down, alternating above and below some value near $0.74$. But writing one more line for every additional step is not a method. The same computation belongs in a loop.
x = 0.3
for k in range(60):
x = np.cos(x)
print(x)
Keeping the history. Overwriting x at every step throws away the information we most want to look at, which is how the iterate approaches its limit. Storing the sequence costs one array.
n = 20
x = np.zeros(n)
x[0] = 0.3
for k in range(n-1):
x[k+1] = np.cos(x[k])
for k in range(n):
print('k = %2d x = %.10f' % (k, x[k]))
Two features of this sequence are worth naming.
The iterates alternate around the solution: one is above, the next below. That is the sign of $g'$ at work. Here $g'(x) = -\sin(x)$ is negative near the root, so $e_{k+1} = g'e_k$ flips the sign of the error at every step.
The approach is steady rather than dramatic: each step removes a fixed fraction of the remaining error. That fraction should be $\lvert g'(r)\rvert = \lvert\sin(r)\rvert$, and the next cell checks it.
r = 0.7390851332151607 # the fixed point, accurate to machine precision
print('r = %.10f' % r)
print('cos(r) = %.10f' % np.cos(r))
print('g\'(r) = %.6f' % (-np.sin(r)))
print('|g\'(r)| = %.6f -> less than 1, so the iteration converges' % abs(np.sin(r)))
print('\n k x_k e_k |e_k| / |e_k-1|')
x = 0.3
e_old = abs(x - r)
for k in range(1, 11):
x = np.cos(x)
e = abs(x - r)
print('%2d %.10f %.3e %.4f' % (k, x, e, e/e_old))
e_old = e
The ratio of successive errors hovers around $0.674$, which is exactly $\lvert g'(r)\rvert$. The theory of Section 1.3.2 is confirmed, and it also tells us the cost: since $0.674^{k}$ must reach $10^{-10}$, about $58$ iterations are needed for full double precision, and that number could have been predicted before running anything.
The cobweb diagram shows the same behaviour geometrically.
plt.figure(figsize = (5, 5))
cobweb(np.cos, 0.3, 12, -0.2, 1.4, 'g(x) = cos(x)')
plt.show()
The path spirals inward toward the crossing point. Spiralling rather than climbing is the geometric form of the sign alternation noticed above, and the fact that it winds inward is the geometric form of $\lvert g'\rvert < 1$: the curve is flatter than the diagonal where they meet.
Does the starting point matter? The convergence result of Section 1.3.2 is local, so it is fair to ask what happens far from the root.
x = 10.0
for k in range(60):
x = np.cos(x)
print(x)
It converges to the same value. This is a happy accident of this particular $g$: no matter where the iteration starts, the first step lands inside $[-1,1]$, because that is the range of the cosine, and from there $\lvert g'\rvert = \lvert\sin(x)\rvert \le \sin(1) \approx 0.84 < 1$ everywhere. The map is a contraction on an interval it cannot escape, so it converges from any starting point. Most fixed-point iterations are not this forgiving.
1.3.4. Example: Computing $\sqrt{2}$¶
Now take $f(x) = x^2 - 2$, whose positive root is $\sqrt{2}$. Rewriting as a fixed-point problem, the most obvious choice is to divide,
$$ x^2 = 2 \quad\Longrightarrow\quad x = \frac{2}{x} = g(x) . $$
This is a correct rewriting: $\sqrt{2}$ certainly satisfies it. Let us iterate.
x = 2.0
for k in range(10):
x = 2/x
print('k = %2d x = %.10f' % (k + 1, x))
The iteration is trapped in a cycle: $2 \to 1 \to 2 \to 1 \to \cdots$, and it never approaches $\sqrt{2} = 1.4142\ldots$ at all. Section 1.3.2 explains why immediately:
$$ g(x) = \frac{2}{x} \qquad\Longrightarrow\qquad g'(x) = -\frac{2}{x^2} \qquad\Longrightarrow\qquad g'(\sqrt{2}) = -\frac{2}{2} = -1 . $$
The magnitude is exactly $1$, the borderline case. The error is neither shrunk nor enlarged, only flipped in sign, so the iteration circles the solution forever without approaching it. The cobweb diagram makes this unmistakable.
x = np.linspace(0.2, 4, 100)
g = lambda x: 2/x
cobweb(g, 2, 10, 0.2, 4, 'g(x) = 2/x')
plt.axis('equal')
plt.ylim([0, 4])
plt.show()
The cobweb closes into a rectangle instead of spiralling in. Nothing is wrong with the algebra; the rewriting is simply a bad one.
A better rewriting. The failure suggests its own cure. When the iterate is too small, $2/x$ overshoots to a value that is too large, and the other way round, so the true answer always lies between $x$ and $2/x$. Taking their average instead of jumping the whole way gives
$$ g(x) = \frac{1}{2}\left(x + \frac{2}{x}\right). $$
This is the same fixed point, since $x = (x + 2/x)/2$ still implies $x^2 = 2$, but a very different iteration.
x = 3.0
for k in range(8):
x = (x + 2/x)/2
print('k = %2d x = %.16f' % (k + 1, x))
print('\nnumpy sqrt(2) = %.16f' % np.sqrt(2))
Six iterations reach machine precision, agreeing with NumPy's own square root to the last stored digit, from a starting point that was not even close. The reason is again in the derivative:
$$ g(x) = \frac{1}{2}\left(x + \frac{2}{x}\right) \qquad\Longrightarrow\qquad g'(x) = \frac{1}{2}\left(1 - \frac{2}{x^2}\right) \qquad\Longrightarrow\qquad g'(\sqrt{2}) = \frac{1}{2}(1-1) = 0 . $$
The derivative vanishes at the root. The linear term in the error recursion disappears, the quadratic term takes over, and the error is squared at every step rather than multiplied by a constant. That is why the number of correct digits roughly doubles per line in the output above: $1, 3, 6, 13$ and then machine precision.
Remark. This is not a lucky guess. Applying Newton's method to $f(x) = x^2-2$ gives $x_{k+1} = x_k - f(x_k)/f'(x_k) = x_k - (x_k^2-2)/(2x_k) = (x_k + 2/x_k)/2$, exactly the averaged iteration. Newton's method is a fixed-point iteration whose $g$ has been constructed so that $g'(r) = 0$, which is precisely why it converges quadratically. The averaging that looked like an ad hoc repair is the systematic construction in disguise.
g_newton = lambda x: (x + 2/x)/2
cobweb(g_newton, 3.0, 6, 0.5, 3.5, 'g(x) = (x + 2/x)/2')
plt.axis('equal')
plt.show()
The staircase descends monotonically onto the crossing point, and the steps shorten so quickly that the last ones are invisible. Compare this with the closed rectangle of $g(x) = 2/x$: the same equation, the same solution, two different rewritings, and completely different behaviour.
1.3.5. Fixed-Point Iteration for Linear Systems¶
The same idea applies to systems of equations. Consider
$$ \begin{align*} 4x_1 - x_2 + x_3 &= 7 \\ 4x_1 - 8x_2 + x_3 &= -21 \\ -2x_1 + x_2 + 5x_3 &= 15 \end{align*} $$
For a system this small the direct solution is immediate.
A = np.array([[4, -1, 1],
[4, -8, 1],
[-2, 1, 5]])
b = np.array([7, -21, 15])
x = np.linalg.solve(A, b)
print(x)
Direct solution is the right tool here, but it costs on the order of $n^3$ operations and it requires storing the whole matrix. For the large sparse systems that come out of discretized PDEs, with $n$ in the millions, iterative methods are preferred, and the simplest of them is the fixed-point iteration of this section.
To build it, solve each equation for the variable on its own diagonal:
$$ \begin{align*} 4x_1 - x_2 + x_3 = 7 \qquad &\Longrightarrow \qquad x_1 = \tfrac{1}{4}x_2 - \tfrac{1}{4}x_3 + \tfrac{7}{4}\\ 4x_1 - 8x_2 + x_3 = -21 \qquad &\Longrightarrow \qquad x_2 = \tfrac{1}{2}x_1 + \tfrac{1}{8}x_3 + \tfrac{21}{8}\\ -2x_1 + x_2 + 5x_3 = 15 \qquad &\Longrightarrow \qquad x_3 = \tfrac{2}{5}x_1 - \tfrac{1}{5}x_2 + 3 \end{align*} $$
Each equation now expresses one unknown in terms of the others, which is exactly the form $\mathbf{x} = G\mathbf{x} + \mathbf{c}$:
$$ \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix} 0 & \frac{1}{4} & -\frac{1}{4}\\ \frac{1}{2} & 0 & \frac{1}{8}\\ \frac{2}{5} & -\frac{1}{5} & 0 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} + \begin{bmatrix} \frac{7}{4}\\ \frac{21}{8}\\ 3 \end{bmatrix} $$
and the iteration is the vector version of $x_{k+1} = g(x_k)$:
$$ \mathbf{x}^{(k+1)} = G\,\mathbf{x}^{(k)} + \mathbf{c} . $$
This particular choice, isolating each variable using its own diagonal entry, is known as the Jacobi method.
G = np.array([[0, 1/4, -1/4],
[1/2, 0, 1/8],
[2/5, -1/5, 0 ]])
c = np.array([7/4, 21/8, 3])
x = np.array([1.0, 1.0, 2.0]) # initial guess
for k in range(20):
x = G.dot(x) + c
print('k = %2d x = %s' % (k + 1, np.round(x, 8)))
The iterates settle onto $(2, 4, 3)$, the same answer the direct solve produced, using nothing but one matrix-vector product per step.
1.3.6. Convergence for Linear Iterations¶
For the linear case the convergence condition can be derived exactly, without the mean value theorem. Unroll the iteration:
$$ \begin{align*} \mathbf{x}^{(k+1)} &= G\mathbf{x}^{(k)} + \mathbf{c}\\ &= G\left(G\mathbf{x}^{(k-1)} + \mathbf{c}\right) + \mathbf{c} = G^{2}\mathbf{x}^{(k-1)} + G\mathbf{c} + \mathbf{c}\\ &\;\;\vdots\\ &= G^{k+1}\mathbf{x}^{(0)} + \left(G^{k} + \cdots + G + I\right)\mathbf{c} . \end{align*} $$
Everything hinges on the powers $G^{k}$. If they shrink to zero, the influence of the initial guess disappears and the second group converges to $(I-G)^{-1}\mathbf{c}$, which is the solution. If they grow, nothing converges.
Writing $G$ in terms of its eigenvalues shows what controls the powers: applying $G$ multiplies the component along each eigenvector by the corresponding eigenvalue $\lambda_i$, so applying it $k$ times multiplies that component by $\lambda_i^{k}$. Every component decays if and only if every eigenvalue has magnitude below one. Defining the spectral radius as the largest of those magnitudes,
$$ \rho(G) = \max_i \lvert \lambda_i \rvert , $$
the condition is
$$ \;\text{the iteration } \quad \mathbf{x}^{(k+1)} = G\mathbf{x}^{(k)} + \mathbf{c} \quad \text{ converges if and only if } \quad \rho(G) < 1\; $$
This is the matrix version of $\lvert g'(r)\rvert<1$, and $\rho(G)$ plays the role of $\lvert g'\rvert$: it is also the factor by which the error shrinks per iteration.
eigenvalues = np.linalg.eigvals(G)
print('eigenvalues of G :', np.round(eigenvalues, 4))
print('spectral radius : %.4f -> less than 1, so it converges' % max(abs(eigenvalues)))
The spectral radius is $0.33$, so each iteration removes about two thirds of the error, which matches the speed observed in Section 1.3.5.
A different rearrangement of the same system. Just as $\sqrt{2}$ had a good rewriting and a bad one, a linear system can be rearranged in more than one way. Using $x = x + f(x)$ or $x = x - f(x)$ row by row instead of dividing by the diagonal gives
$$ \begin{align*} x_1 &= -3x_1 + x_2 - x_3 + 7\\ x_2 &= 4x_1 - 7x_2 + x_3 + 21\\ x_3 &= 2x_1 - x_2 - 4x_3 + 15 \end{align*} $$
Substituting $(2,4,3)$ confirms that this has the same solution as before. Now iterate it.
G2 = np.array([[-3, 1, -1],
[ 4, -7, 1],
[ 2, -1, -4]])
c2 = np.array([7, 21, 15])
x = np.array([1.0, 2.0, 2.0])
for k in range(10):
x = G2.dot(x) + c2
print('k = %2d x = %s' % (k + 1, np.round(x, 2)))
The solution is unchanged, and the iteration is useless: within ten steps the iterates have passed $10^{8}$. The spectral radius says why.
eigenvalues2 = np.linalg.eigvals(G2)
print('eigenvalues of G2 :', np.round(eigenvalues2, 4))
print('spectral radius : %.4f -> greater than 1, so it diverges' % max(abs(eigenvalues2)))
Every eigenvalue lies outside the unit circle, the largest with magnitude $7.16$, so each iteration multiplies the error by roughly seven. After ten steps that is a factor of $7^{10}\approx 3\times10^{8}$, which is what the output shows.
1.4.1. The Idea: Follow the Tangent¶
Section 1.3 solved $f(x)=0$ by rewriting it as $x=g(x)$ and iterating, and found that everything depended on a choice we made by hand. A good $g$ converged quickly, a bad one cycled forever, and there was no systematic recipe for producing a good one.
Newton's method is that recipe. Instead of guessing a rearrangement, it builds one from the function itself, using the only information available at a point: the value there and the slope there.
Suppose the current guess is $x_k$. Expand $f$ about that point and keep the linear term,
$$ f(x) \approx f(x_k) + f'(x_k)\,(x - x_k) . $$
This straight line is the tangent to $f$ at $x_k$. We cannot solve $f(x)=0$, but we can certainly solve the tangent equal to zero:
$$ f(x_k) + f'(x_k)\,(x_{k+1} - x_k) = 0 \qquad\Longrightarrow\qquad \;x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}\; $$
The geometric statement. Stand on the curve at $x_k$, slide down the tangent line until it crosses the axis, and take that crossing as the next guess. Repeat. Where the curve is well approximated by its tangent, which is to say near a simple root, the crossing of the tangent is very close to the crossing of the curve, and the method converges extremely fast.
The price. The formula needs $f'$, which the fixed-point iteration of Section 1.3 did not. It must be derived by hand, obtained by automatic differentiation, or approximated. It must also not vanish: the method divides by it.
1.4.2. Newton's Method Is a Fixed-Point Iteration¶
Newton's method is not a new kind of algorithm. Reading the update as $x_{k+1}=g(x_k)$ identifies its map,
$$ g(x) = x - \frac{f(x)}{f'(x)} , $$
so everything proved in Section 1.3.2 applies. In particular the speed of convergence is decided by $g'$ at the root, and that derivative is worth computing. Differentiating with the quotient rule,
$$ g'(x) = 1 - \frac{f'(x)f'(x) - f(x)f''(x)}{\left[f'(x)\right]^2} = \frac{f(x)\,f''(x)}{\left[f'(x)\right]^2} . $$
Now evaluate at a simple root $r$, meaning one where $f(r)=0$ and $f'(r)\neq0$. The numerator contains the factor $f(r)$, which is zero, and the denominator is not zero. Therefore
$$ g'(r) = 0\; $$
This single line explains everything about the method. Section 1.3.2 observed that $g'(r)=0$ is the exceptional case in which the linear error term disappears and convergence becomes faster than linear. Newton's method is the construction that produces that case on purpose, for any $f$. The averaging trick that repaired the $\sqrt{2}$ iteration in Section 1.3.4 was not an isolated piece of cleverness; it was this formula applied to $f(x)=x^2-2$.
Notice also what the derivation assumed. If $f'(r)=0$, the root is not simple, the cancellation above fails, and the argument collapses. At such a root the method still converges, but only linearly.
1.4.3. Quadratic Convergence¶
To see how fast "faster than linear" is, expand the error one term further than Section 1.3.2 did. Write $e_k = x_k - r$ and expand $f$ and $f'$ about the root:
$$ f(x_k) = f'(r)\,e_k + \tfrac{1}{2}f''(r)\,e_k^2 + \cdots, \qquad f'(x_k) = f'(r) + f''(r)\,e_k + \cdots . $$
Substituting into $e_{k+1} = e_k - f(x_k)/f'(x_k)$ and keeping terms through $e_k^2$ gives
$$ e_{k+1} \approx \frac{f''(r)}{2f'(r)}\;e_k^{2} $$
The error is squared at every step, not multiplied by a constant. This is quadratic convergence, and its practical meaning is that the number of correct digits doubles per iteration: a guess good to two digits becomes good to four, then eight, then sixteen, at which point double precision is exhausted. Three or four iterations from a reasonable starting point is typical, and it is why Newton's method is the default for nonlinear problems throughout computational engineering.
The constant in front is not negligible but it is rarely decisive, because squaring the error overwhelms it after a step or two. It does say something useful, though: convergence is fastest where the curve is straightest near the root, since $f''(r)$ small means the tangent line is a good model over a wide neighbourhood.
The test problem. Take the transcendental equation
$$ f(x) = x\,e^{x} - 2 = 0 . $$
It cannot be solved by algebra: the unknown appears both as a factor and inside an exponent, and no rearrangement separates the two. Equations of this shape are common in engineering, arising whenever a linear effect multiplies an exponential one, as in thermal, chemical-kinetic and semiconductor problems.
The derivatives follow from the product rule,
$$ f'(x) = e^{x} + x\,e^{x} = (1+x)\,e^{x}, \qquad f''(x) = (2+x)\,e^{x} , $$
and both are positive for $x>0$, so there $f$ is increasing and convex. Since $f(0)=-2<0$ and $f(1)=e-2>0$, there is exactly one positive root, lying between $0$ and $1$. Its exact value is the Lambert $W$ function evaluated at $2$, $r = W(2) = 0.8526055020137254$, which we use only as a reference for measuring error.
The error law predicts the constant
$$ \frac{f''(r)}{2f'(r)} = \frac{(2+r)e^{r}}{2(1+r)e^{r}} = \frac{2+r}{2(1+r)} = 0.7699 , $$
the exponential cancelling. The next cell starts deliberately far away, at $x_0 = 2$.
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return x*np.exp(x) - 2
def dfdx(x):
return (1 + x)*np.exp(x)
r = 0.8526055020137254 # the root, W(2), accurate to machine precision
x = 2.0
e_old = abs(x - r)
print(' k x_k e_k e_k / e_(k-1)^2')
for k in range(1, 9):
x = x - f(x)/dfdx(x)
e = abs(x - r)
if e > 0:
print('%2d %.16f %.3e %.4f' % (k, x, e, e/e_old**2))
e_old = e
else:
print('%2d %.16f %.3e (exact)' % (k, x, e))
print('\npredicted constant f\'\'(r) / (2 f\'(r)) = %.4f' % ((2 + r)/(2*(1 + r))))
Two phases are visible, and the separation between them is the lesson.
The first three steps are not fast. Starting at $x_0=2$ the error falls only from $1.15$ to $5.7\times10^{-1}$, $1.8\times10^{-1}$, $2.3\times10^{-2}$, roughly a factor of three per step, because the tangent at $x=2$ is a poor model of a curve that is bending sharply there. This is the approach phase, and no theory of this section governs it.
From the fourth step on the behaviour changes completely: $4.0\times10^{-4}$, $1.2\times10^{-7}$, $1.2\times10^{-14}$, then exact. Two correct digits become four, then seven, then fourteen, then all of them. Once the iterate is close enough for the quadratic law to hold, four steps exhaust double precision.
The ratio column confirms the law numerically. It drifts upward through the approach phase and then settles on $0.7697$, matching the predicted $(2+r)/\bigl(2(1+r)\bigr)=0.7699$ to three digits. That constant was computed from the analytic derivatives alone, before any iteration was run.
Comparison with a fixed-point iteration. The same equation can be rearranged in the manner of Section 1.3. Dividing $x e^{x}=2$ by $e^{x}$ gives
$$ x = 2e^{-x} \;\equiv\; g(x), \qquad g'(x) = -2e^{-x} . $$
At the root, $2e^{-r}=r$ by the equation itself, so $g'(r) = -r = -0.8526$. Its magnitude is less than one, so by the criterion of Section 1.3.2 the iteration converges, and the negative sign means it alternates about the root. But $0.85$ is uncomfortably close to $1$: each step removes only $15\%$ of the error, so roughly fifteen iterations are needed per decimal digit.
# fixed-point iteration x <- 2 exp(-x)
x = 2.0
fixed_point_errors = []
for k in range(10):
x = 2*np.exp(-x)
fixed_point_errors.append(abs(x - r))
# Newton's method on f(x) = x exp(x) - 2
x = 2.0
newton_errors = []
for k in range(10):
x = x - f(x)/dfdx(x)
newton_errors.append(abs(x - r))
floor = 1e-17 # so that exact zeros can be drawn on a log axis
plt.figure(figsize = (6, 4))
plt.semilogy(np.maximum(fixed_point_errors, floor), 'b', linewidth = 2, label = 'fixed point : x = 2 exp(-x)')
plt.semilogy(np.maximum(newton_errors, floor), 'rs-', label = 'Newton : f(x) = x exp(x) - 2')
plt.xlabel('iteration k'); plt.ylabel('|error|'); plt.legend(); plt.show()
On a logarithmic axis, linear convergence is a straight line and quadratic convergence is a curve that bends downward ever more steeply. The fixed-point iteration descends at a constant slope and needs $221$ steps to reach machine precision; Newton's method falls off the bottom of the plot in seven, so far to the left that its curve is pressed against the axis. The two methods are solving the same equation from the same starting point.
1.4.4. The Geometric Picture¶
The algebra is short enough that the geometry is easy to lose sight of. The next figure draws the first three tangent lines explicitly, starting again from $x_0 = 2$.
xs = np.linspace(0.4, 2.2, 200)
plt.figure(figsize = (6, 4))
plt.plot(xs, f(xs), 'k', linewidth = 2, label = r'$f(x) = x e^x - 2$')
plt.axhline(y = 0, color = 'k', linewidth = 1)
x = 2.0
for k in range(3):
slope = dfdx(x)
x_new = x - f(x)/slope
seg = np.linspace(x_new - 0.15, x + 0.1, 20)
plt.plot(seg, f(x) + slope*(seg - x), 'r', alpha = 0.7)
plt.plot([x, x], [0, f(x)], 'k--', alpha = 0.5, linewidth = 1)
plt.plot(x, f(x), 'ro', markersize = 5)
plt.plot(x_new, 0, 'r')
plt.text(x + 0.03, -1.0, '$x_%s$' % k)
x = x_new
plt.plot(r, 0, 'ko', markersize = 5)
plt.xlim([0.4, 2.2]); plt.ylim([-3, 14])
plt.xlabel('x'); plt.ylabel('f(x)'); plt.legend(); plt.show()
Read the figure as a sequence. Each circle is the current iterate lifted onto the curve, each red line is the tangent there, and each cross is where that tangent meets the axis, which becomes the next iterate.
The picture explains both phases of the table. Far to the right the curve is strongly convex, the tangent lies well below it, and the crossing at $x_1 = 1.424$ falls far short of the root. The second tangent, taken on a gentler part of the curve, lands at $x_2 = 1.035$, much closer. By the third the curve is nearly straight over the remaining distance, the tangent is almost indistinguishable from it, and the crossing at $x_3 = 0.876$ already sits on the star.
Two features of the geometry are worth naming. Because $f$ is convex and increasing here, every tangent line lies below the curve, so every crossing falls to the right of the root and the iterates approach from one side only. And the horizontal distances from the crosses to the star, $0.57$, $0.18$, $0.023$, are exactly the errors of the table: the figure and the table are two views of the same three steps.
1.4.5. Newton's Method for Systems¶
Nothing in the derivation was one-dimensional. Given $n$ equations in $n$ unknowns,
$$ \mathbf{F}(\mathbf{x}) = \mathbf{0}, \qquad \mathbf{F} : \mathbb{R}^n \to \mathbb{R}^n , $$
expand about the current iterate and keep the linear term. The role of the derivative is taken by the Jacobian matrix, the matrix of all first partial derivatives,
$$ J_{ij}(\mathbf{x}) = \frac{\partial F_i}{\partial x_j} , \qquad \mathbf{F}(\mathbf{x}) \approx \mathbf{F}(\mathbf{x}_k) + J(\mathbf{x}_k)\,(\mathbf{x}-\mathbf{x}_k) , $$
and setting the linear model to zero gives the multivariable update
$$ J(\mathbf{x}_k)\,\Delta\mathbf{x} = -\mathbf{F}(\mathbf{x}_k), \qquad \mathbf{x}_{k+1} = \mathbf{x}_k + \Delta\mathbf{x} $$
Two points deserve emphasis. Division by $f'$ has become the solution of a linear system, so every Newton step for a system costs one matrix factorization, which is where the linear solvers of the numerical analysis material earn their place. And one should never form $J^{-1}$: solving $J\,\Delta\mathbf{x}=-\mathbf{F}$ is both cheaper and more accurate than computing an inverse and multiplying by it.
Test problem. Find where a circle and a parabola intersect,
$$ \begin{align*} x^2 + y^2 &= 4\\ y &= x^2 \end{align*} \qquad\Longrightarrow\qquad \mathbf{F}(x,y)=\begin{bmatrix} x^2+y^2-4 \\ y-x^2 \end{bmatrix}, \qquad J=\begin{bmatrix} 2x & 2y \\ -2x & 1 \end{bmatrix}. $$
def F(v):
x, y = v
return np.array([x**2 + y**2 - 4,
y - x**2])
def Jacobian(v):
x, y = v
return np.array([[ 2*x, 2*y],
[-2*x, 1.0]])
v = np.array([1.0, 1.0]) # initial guess
print(' k x y |F(v)|')
for k in range(5):
delta = np.linalg.solve(Jacobian(v), -F(v))
v = v + delta
print('%2d %.10f %.10f %.3e' % (k + 1, v[0], v[1], np.linalg.norm(F(v))))
y_exact = (-1 + np.sqrt(17))/2
print('\nexact : x = %.10f y = %.10f' % (np.sqrt(y_exact), y_exact))
The residual falls as $5.7\times10^{-1}$, $1.8\times10^{-2}$, $2.3\times10^{-5}$, $4.9\times10^{-11}$, and then round-off: the same doubling of correct digits as in one dimension. Eliminating $y$ by hand gives $x^2+x^4=4$, so the exact intersection has $y=(\sqrt{17}-1)/2$ and $x=\sqrt{y}$, which the last line confirms.
theta = np.linspace(0, 2*np.pi, 200)
xs = np.linspace(-2.2, 2.2, 200)
plt.figure(figsize = (5, 5))
plt.plot(2*np.cos(theta), 2*np.sin(theta), 'k', linewidth = 1, alpha = 0.6, label = r'$x^2 + y^2 = 4$')
plt.plot(xs, xs**2, 'k', linewidth = 1, alpha = 0.6, label = r'$y = x^2$')
v = np.array([1.0, 0.5])
path = [v]
for k in range(2):
v = v + np.linalg.solve(Jacobian(v), -F(v))
path.append(v)
path = np.array(path)
plt.plot(path[:, 0], path[:, 1], 'ro', markersize = 3, label = 'Newton iterates')
plt.xlim([-2.5, 2.5]); plt.ylim([-2.5, 2.5])
plt.grid(alpha = 0.3)
plt.xlabel('x'); plt.ylabel('y'); plt.legend(); plt.show()
The iterates start off both curves and move onto their intersection within two visible steps. Which of the two intersections is found depends entirely on the starting guess; beginning at $(-1,1)$ would have produced the mirror solution. Newton's method finds a root, never the roots, and in a system with many solutions the initial guess is what selects among them.
1.4.6. Bisection, Fixed Point, and Newton¶
Newton's method answers the question left open at the end of Section 1.3: how to choose the fixed-point map $g$ rather than guess it. The answer, $g = x - f/f'$, is constructed so that $g'(r)=0$ at every simple root, which is exactly the condition that makes convergence quadratic instead of linear.
What is gained is speed, and the margin is not small. On the same equation, plain fixed-point iteration needed about sixty steps where Newton needed four, and in the system of Section 1.4.5 the residual went from $10^{-1}$ to round-off in four steps.
What is given up is robustness. Fixed-point iteration with $\lvert g'\rvert<1$ converges from anywhere in the region where that bound holds; Newton's method carries no such guarantee: its convergence is only local, a step landing where $f'$ is near zero can be thrown far across the domain, and at a multiple root the rate drops back to linear. It also requires a derivative, or a Jacobian for a system, and a linear solve at every step.
The practical conclusion is the one every numerical library implements: use a slow and safe method to get close, then switch to Newton to finish. Bisection to locate, Newton to converge.
1.5.1. Derivation and Implementation¶
Newton's method needs $f'$, which may be unavailable, expensive, or tedious to derive. The secant method replaces the tangent line by the line through the last two points $(x_{k-1}, f(x_{k-1}))$ and $(x_k, f(x_k))$, i.e. it approximates $f'(x_k)$ by a difference quotient:
$$ \begin{align*} f'(x_k) &\approx \frac{f(x_k)-f(x_{k-1})}{x_k - x_{k-1}}\\\\ x_{k+1} &= x_k - \frac{f(x_k)\,(x_k - x_{k-1})}{f(x_k) - f(x_{k-1})}, \qquad k = 2, 3, \dots \end{align*} $$
Two starting values $x_1, x_2$ are needed, and each iteration costs only one new function evaluation, since $f(x_{k-1})$ is remembered from the previous step.
Because the secant slope is only an approximation of the tangent slope, the convergence is not quite quadratic; the error obeys $|\epsilon_{k+1}| \approx c\,|\epsilon_k|^{\alpha}$ with $\alpha = (1 + \sqrt 5)/2 \approx 1.618$, the golden ratio. This is called superlinear convergence: slower than Newton per iteration, but each iteration is half the work, since Newton needs both $f$ and $f'$.
x1 = 1.0
x2 = 0.5
xs = [x1, x2]
# take secant steps until the residual is negligible
while abs(f(x2)) > 1e-12:
x3 = x2 - f(x2)*(x2 - x1)/(f(x2) - f(x1))
x1 = x2
x2 = x3
xs.append(x2)
xs = np.array(xs)
errs = np.abs(xs - r)
print(' k x_k |error| log|e_k+1|/log|e_k|')
for k in range(len(xs)):
if k + 1 < len(xs):
ratio = np.log(errs[k+1])/np.log(errs[k])
print('%3d %.12f %9.2e %6.3f' % (k + 1, xs[k], errs[k], ratio))
else:
print('%3d %.12f %9.2e' % (k + 1, xs[k], errs[k]))
2.1. Why Numerical Derivatives?¶
2.1.1. Motivation¶
Differentiation is taught as a symbolic operation: given a formula for $f$, the rules of calculus produce a formula for $f'$. On a computer that route is frequently unavailable.
- The function is known only at sampled points — a sensor record, a table of measurements, the output of an experiment. There is no formula to differentiate.
- The function is the output of a simulation or a black box — a finite-element solver, a legacy subroutine — that returns a number for each input and nothing else.
- The derivative is needed at many points of a grid, as the building block of a numerical method for a differential equation. Discretizing $\partial u/\partial t = \alpha\, \partial^2 u/\partial x^2$ means replacing both derivatives by differences of grid values; this is the origin of the name finite difference method.
In each case we need a rule that turns values of $f$ into an approximation of $f'$, together with an honest statement of how wrong that approximation is.
2.1.2. The Limit That Cannot Be Taken¶
The definition of the derivative,
$$ f'(x) = \lim_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x)}{\Delta x} , $$
already suggests the approximation: stop before the limit and use a small but finite $\Delta x$. That single decision — the finite in finite difference — is the source of everything that follows.
Two questions arise immediately, and they pull in opposite directions. How large is the error committed by stopping at a finite $\Delta x$? It shrinks as $\Delta x$ shrinks, so we want $\Delta x$ small. But the numerator $f(x + \Delta x) - f(x)$ is a difference of two nearly equal numbers, and in floating-point arithmetic such a difference loses significant digits; dividing the damaged numerator by a tiny $\Delta x$ magnifies the damage. So we also want $\Delta x$ not small. The sections that follow study the first of these two errors.
Throughout we use the test function
$$ f(x) = \sin x , \qquad f'(x) = \cos x , \qquad f''(x) = -\sin x , $$
whose derivatives are known exactly, so that every approximation can be compared with the truth. All derivatives of $\sin$ are bounded by $1$, which will make the error estimates especially clean.
2.2.1. Derivation from Taylor Series¶
Everything in this chapter comes from two expansions, one step to the right and one step to the left:
$$ (*)\qquad f(x + \Delta x) = f(x) + \Delta x\, f'(x) + \frac{\Delta x^2}{2!}f''(x) + \frac{\Delta x^3}{3!}f'''(x) + \cdots $$
$$ (**)\qquad f(x - \Delta x) = f(x) - \Delta x\, f'(x) + \frac{\Delta x^2}{2!}f''(x) - \frac{\Delta x^3}{3!}f'''(x) + \cdots $$
Solving $(*)$ for $f'(x)$ and dividing by $\Delta x$ gives the forward difference,
$$ \frac{f(x + \Delta x) - f(x)}{\Delta x} = f'(x) + \frac{\Delta x}{2!}f''(x) + \frac{\Delta x^2}{3!}f'''(x) + \cdots $$
The expression tells us both that the left-hand side approximates $f'(x)$ and by how much it misses:
$$ D_{+}f(x) = \frac{f(x + \Delta x) - f(x)}{\Delta x}, \qquad D_{+}f(x) - f'(x) = \frac{\Delta x}{2}f''(x) + O(\Delta x^2) . $$
The error is proportional to the first power of $\Delta x$; the formula is said to be first-order accurate, or of order $O(\Delta x)$. Halving $\Delta x$ halves the error.
The same treatment of $(**)$ gives the backward difference
$$ D_{-}f(x) = \frac{f(x) - f(x - \Delta x)}{\Delta x}, \qquad D_{-}f(x) - f'(x) = -\frac{\Delta x}{2}f''(x) + O(\Delta x^2) , $$
also first-order accurate. The two leading error terms are equal in size and opposite in sign — a fact that Section 2.3 exploits.
Remark (geometry). The forward difference is the slope of the line joining $(x, f(x))$ to the point one step to the right; the backward difference uses the point one step to the left. The tangent line at $x$ lies between them, which is exactly what the opposite signs of the two error terms say.
2.2.2. The Two Slopes at a Point¶
The three formulas are evaluated below at $x_0 = 1$ with a deliberately large step $\Delta x = 0.25$, so that the differences between them are easy to see. The figure draws the three sample points $x - \Delta x$, $x$ and $x + \Delta x$ together with the central-difference line through the outer two.
x0 = 1
dx = .25
# Forward Difference Approximation
dfdxF = (np.sin(x0 + dx) - np.sin(x0)) / dx
# Backward Difference Approximation
dfdxB = (np.sin(x0) - np.sin(x0 - dx)) / dx
dfdxC = (np.sin(x0 + dx) - np.sin(x0 - dx)) / (2*dx)
print('exact %8.5f' % np.cos(x0))
print('forward %8.5f error %8.5f' % (dfdxF, dfdxF - np.cos(x0)))
print('backward %8.5f error %8.5f' % (dfdxB, dfdxB - np.cos(x0)))
# plotting commands
xs = np.linspace(0.3, 1.7, 200)
plt.figure(figsize = (6, 4))
plt.plot(xs, np.sin(xs), 'k', linewidth = 1)
plt.plot(xs, np.sin(x0) + np.cos(x0)*(xs - x0), 'k--', linewidth = 1.2) # exact tangent
plt.plot(xs, np.sin(x0) + dfdxF*(xs - x0), 'b', linewidth = 1.2) # Forward Difference
plt.plot(xs, np.sin(x0) + dfdxB*(xs - x0), 'g', linewidth = 1.2) # Backward Difference
# plt.plot(xs, np.sin(x0) + dfdxC*(xs - x0), 'm', linewidth=1.2) # Central Difference
plt.plot([x0 - dx, x0, x0 + dx], np.sin([x0 - dx, x0, x0 + dx]), 'ro', markersize = 4)
plt.plot([x0 - dx, x0 - dx], [0, np.sin(x0 - dx)], 'k--', alpha = 0.5, linewidth = 1)
plt.plot([x0, x0], [0, np.sin(x0)], 'k--', alpha = 0.5, linewidth = 1)
plt.plot([x0 + dx, x0 + dx], [0, np.sin(x0 + dx)], 'k--', alpha = 0.5, linewidth = 1)
plt.text(x0 + 0.03, 0.52, '$x$')
plt.text(x0 - dx + 0.03, 0.52, r'$x - \Delta x$')
plt.text(x0 + dx + 0.03, 0.52, r'$x + \Delta x$')
# plt.legend(['Function', 'Exact Tangent', 'Forward Diff', 'Backward Diff'])
plt.ylim(0.5, 1.15)
plt.show()
The forward slope is too small by $0.110$ and the backward slope too large by $0.099$ — opposite in sign and both close to the predicted magnitude $\tfrac{\Delta x}{2}|f''(x_0)| = 0.105$. The true tangent therefore lies between the two secant slopes, as the remark predicted, and the central difference — the line drawn in the figure — is their average.
2.3.1. Derivation¶
The forward and backward differences have leading errors $+\tfrac{\Delta x}{2}f''$ and $-\tfrac{\Delta x}{2}f''$. Averaging the two formulas must therefore cancel that term. Equivalently, and more systematically, subtract $(**)$ from $(*)$. Every term of even order is identical in the two expansions and cancels:
$$ f(x + \Delta x) - f(x - \Delta x) = 2\Delta x\, f'(x) + \frac{2\Delta x^3}{3!}f'''(x) + O(\Delta x^5) $$
Dividing by $2\Delta x$ gives the central difference
$$ D_{0}f(x) = \frac{f(x + \Delta x) - f(x - \Delta x)}{2\Delta x}, \qquad D_{0}f(x) - f'(x) = \frac{\Delta x^2}{3!}f'''(x) + O(\Delta x^4) . $$
The formula is second-order accurate: halving $\Delta x$ divides the error by four. It costs the same two function evaluations as the forward difference, so the extra accuracy is free — the only price is that it needs a point on each side, which rules it out in two common situations: when $f'$ must be computed in real time, where the future value $f(x + \Delta x)$ is not yet available, and at the two ends of a finite data record (Section 2.5).
Intuition (why symmetry pays). The cancellation was not an accident of algebra. A symmetric formula cannot distinguish left from right, so its error cannot contain any term that changes sign when $\Delta x \to -\Delta x$ — that is, no odd power of $\Delta x$. Symmetry alone forces the error to start at $\Delta x^2$. The same argument explains why every symmetric stencil below has an even order of accuracy.
2.3.2. Comparing the Three Approximations¶
The tangent picture of Section 2.2.2 was drawn at a single point. A more revealing experiment applies the three formulas at every point of an interval and plots the results against the exact derivative. Take $f = \sin$ on $[-2, 4]$ with a coarse $\Delta t = 0.2$, so that the differences between the schemes are visible to the eye.
dt = .2
t = np.linspace(-2, 4, 60)
f = np.sin(t)
# Exact Derivative
dfdt = np.cos(t)
# plotting commands
plt.plot(t, f, 'k--', linewidth = 1.2)
plt.grid(True)
plt.plot(t, dfdt, 'k', linewidth = 3)
plt.legend(['Function', 'Exact Derivative'], fontsize = 14)
plt.xlim(-2, 4)
plt.ylim(-1.5, 1.5)
plt.show()
# Forward Difference Approximation
dfdtF = (np.sin(t + dt) - np.sin(t)) / dt
# Backward Difference Approximation
dfdtB = (np.sin(t) - np.sin(t - dt)) / dt
# Central Difference Approximation
dfdtC = (np.sin(t + dt) - np.sin(t - dt)) / (2 * dt)
plt.figure(figsize = (20, 4))
plt.subplot(1, 2, 1)
plt.plot(t, dfdt, 'k', linewidth = 3)
plt.plot(t, dfdtF, 'b', linewidth = 1.2) # Forward Difference
plt.plot(t, dfdtB, 'g', linewidth = 1.2) # Backward Difference
plt.plot(t, dfdtC, 'r', linewidth = 1.2) # Central Difference
plt.legend(['Exact Derivative', 'Forward Diff', 'Backward Diff', 'Central Diff'],
fontsize = 14, framealpha = 1.0) # get rid of legend transparency
plt.grid(True)
# now make a zoomed-in plot
plt.subplot(1, 2, 2)
plt.plot(t, dfdt, 'k', linewidth = 3)
plt.plot(t, dfdtF, 'b', linewidth = 1.2) # Forward Difference
plt.plot(t, dfdtB, 'g', linewidth = 1.2) # Backward Difference
plt.plot(t, dfdtC, 'r', linewidth = 1.2) # Central Difference
plt.grid(True)
plt.xlim(0, 2)
plt.show()
# how large is the error of each one, over the whole interval?
print('max error, dt = %g' % dt)
print(' forward %8.5f' % np.max(np.abs(dfdtF - dfdt)))
print(' backward %8.5f' % np.max(np.abs(dfdtB - dfdt)))
print(' central %8.5f' % np.max(np.abs(dfdtC - dfdt)))
# is the forward curve the exact derivative shifted by half a step?
print('\nmax | forward - cos(t + dt/2) | %8.5f' % np.max(np.abs(dfdtF - np.cos(t + dt/2))))
print('max | backward - cos(t - dt/2) | %8.5f' % np.max(np.abs(dfdtB - np.cos(t - dt/2))))
# and how much amplitude does the central difference lose?
print('\namplitude of central / exact %8.5f' % (np.max(np.abs(dfdtC)) / np.max(np.abs(dfdt))))
print('predicted 1 - dt^2/6 %8.5f' % (1 - dt**2/6))
The zoomed panel shows what the error formulas predicted. The forward curve lies to the left of the exact derivative and the backward curve to the right, each by about half a step, and the printed check confirms it: the forward curve agrees with $\cos(t + \Delta t/2)$ sixty times better than with $\cos t$ itself. A forward difference is really the derivative evaluated near the midpoint $t + \Delta t/2$, not at $t$. The central curve is the average of the two, so the shifts cancel; it sits on top of the exact derivative and is hard to see at all, with an error fifteen times smaller than the other two.
Note also that the central curve has a slightly smaller amplitude than $\cos t$, and the measured ratio $0.9934$ is the factor $1 - \Delta t^2/6 = 0.9933$ that the error term predicts. A difference quotient over a finite step averages the function's variation across that step, and averaging always shaves the peaks. This is why finite differences act as mild low-pass filters on the data they are applied to.
2.4.1. The Second Derivative¶
One route to $f''$ is to apply a difference formula twice: $f''$ is the derivative of $f'$, so
$$ f''(x) \approx \frac{f'(x + \Delta x) - f'(x)}{\Delta x} , $$
with each $f'$ replaced in turn by a difference of $f$ values. Any of the three schemes can be used, but it is important to be consistent and use the same scheme for both terms, or the orders of accuracy will not combine as expected.
The direct route is shorter. The central difference came from subtracting the two Taylor expansions; adding them cancels the odd terms instead:
$$ f(x + \Delta x) + f(x - \Delta x) = 2f(x) + \Delta x^2 f''(x) + \frac{\Delta x^4}{4!}f''''(x) + O(\Delta x^6) $$
Solving for $f''(x)$ gives the standard three-point formula
$$ \frac{f(x + \Delta x) - 2f(x) + f(x - \Delta x)}{\Delta x^2} = f''(x) + O(\Delta x^2) , $$
second-order accurate, and exactly what the two-stage route produces when central differences are used at both stages. Written on a grid with $f_k = f(x_k)$ and uniform spacing $\Delta x$, the two workhorse formulas of the finite difference method are
$$ f'(x_k) \approx \frac{f_{k+1} - f_{k-1}}{2\Delta x}, \qquad f''(x_k) \approx \frac{f_{k+1} - 2f_k + f_{k-1}}{\Delta x^2} . $$
The coefficient pattern $(1, -2, 1)/\Delta x^2$ is the discrete Laplacian; it reappears in every finite difference treatment of a diffusion equation, a beam equation or a Poisson problem.
2.4.2. Higher Order by Combining Stencils¶
Accuracy can be raised systematically by using more points. The central difference with step $\Delta x$ and the same formula with step $2\Delta x$ have known error expansions,
$$ D_0(\Delta x) = f' + \frac{\Delta x^2}{6}f''' + O(\Delta x^4), \qquad D_0(2\Delta x) = f' + \frac{4\Delta x^2}{6}f''' + O(\Delta x^4) . $$
The combination $\tfrac{1}{3}\bigl(4 D_0(\Delta x) - D_0(2\Delta x)\bigr)$ is built precisely so that the $\Delta x^2$ terms cancel, leaving
$$ f'(x) \approx \frac{-f(x + 2\Delta x) + 8f(x + \Delta x) - 8f(x - \Delta x) + f(x - 2\Delta x)}{12\,\Delta x} = f'(x) + O(\Delta x^4) . $$
This is the fourth-order central difference. Its five-point stencil uses four function evaluations instead of two, and by symmetry the next error term is $\Delta x^4$, not $\Delta x^3$. The same bookkeeping — write the Taylor expansion of each sample in the stencil, then choose the coefficients so that as many low-order terms as possible cancel — produces a finite difference formula for any derivative to any order of accuracy, and is how the coefficient tables in reference books are generated.
dxs = np.logspace(-1, -5, 17) # step sizes, from coarse to fine
# First Derivative, 3-point stencil
errC = np.abs((np.sin(x0 + dxs) - np.sin(x0 - dxs)) / (2 * dxs) - np.cos(x0))
# Second Derivative, 3-point stencil
err2 = np.abs((np.sin(x0 + dxs) - 2*np.sin(x0) + np.sin(x0 - dxs)) / dxs**2 + np.sin(x0))
# First Derivative, 5-point stencil
dfdx4 = (-np.sin(x0 + 2*dxs) + 8*np.sin(x0 + dxs)
- 8*np.sin(x0 - dxs) + np.sin(x0 - 2*dxs)) / (12 * dxs)
err4 = np.abs(dfdx4 - np.cos(x0))
m = dxs > 1e-2 # fit where truncation error dominates
print('second derivative, 3-point stencil : slope %5.3f' % np.polyfit(np.log10(dxs[m]), np.log10(err2[m]), 1)[0])
print('first derivative, 5-point stencil : slope %5.3f' % np.polyfit(np.log10(dxs[m]), np.log10(err4[m]), 1)[0])
plt.loglog(dxs, errC, 'r^-', linewidth = 1.2)
plt.loglog(dxs, err4, 'md-', linewidth = 1.2)
plt.loglog(dxs, err2, 'cv-', linewidth = 1.2)
plt.legend(["$f'$, 3-point, $O(\\Delta x^2)$", "$f'$, 5-point, $O(\\Delta x^4)$",
"$f''$, 3-point, $O(\\Delta x^2)$"], fontsize = 12)
plt.grid(True, which = 'both')
plt.xlabel('$\\Delta x$', fontsize = 12)
plt.ylabel('absolute error', fontsize = 12)
plt.show()
The measured slopes are $2.00$ and $4.00$. Two curves stop descending before the right edge of the plot: the fourth-order one bottoms out near $\Delta x \approx 5\times10^{-4}$, simply because it has already run out of digits to gain, and the second-derivative one near $\Delta x \approx 5\times10^{-5}$, because dividing by $\Delta x^2$ instead of $\Delta x$ amplifies roundoff far more aggressively.
2.5.1. Interior Points and the Two Ends¶
In practice a function is rarely differentiated at one point. It is sampled on a grid
$$ x_1, x_2, \dots, x_n , \qquad x_k - x_{k-1} = \Delta x , \qquad f_k = f(x_k) , $$
and the whole array of derivatives is wanted at once. Here $\Delta x$ is no longer a free parameter — it is fixed by how densely the data were sampled.
The recipe follows directly from Sections 2.2 and 2.3. Use the central difference on the interior points $x_2, \dots, x_{n-1}$, where a neighbour exists on each side. At $x_1$ there is no point to the left, so use a forward difference; at $x_n$ there is no point to the right, so use a backward difference.
The sampled function is plotted first, with the sample points marked, as a reminder that the grid values are all the information available.
n = 30
x = np.linspace(0.1, 3, n)
f = np.sin(x)
plt.plot(x, f, 'k')
plt.plot(x, f, 'rx', linewidth = 2)
plt.grid(True)
plt.show()
print('grid spacing dx = %.4f' % (x[1] - x[0]))
The differencing itself is a three-line loop. Written this way it follows the recipe literally, one grid point at a time.
$$\frac{f(x_2) - f(x_1)}{\Delta x}$$
$$\frac{f(x_{k+1}) - f(x_{k-1})}{2 \Delta x}$$
$$\frac{f(x_n) - f(x_{n-1})}{\Delta x}$$
dx = x[1] - x[0]
dfdx = np.zeros(n)
# Perform the finite differencing
dfdx[0] = (f[1] - f[0]) / (x[1] - x[0]) # forward diff at f(x_1)
for i in range(1, n - 1):
dfdx[i] = (f[i + 1] - f[i - 1]) / (x[i + 1] - x[i - 1]) # central in between
dfdx[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2]) # backward diff at f(x_n)
plt.plot(x, np.cos(x), 'k', label = 'true derivative')
plt.plot(x, dfdx, 'rx', label = 'computed derivative')
plt.legend(fontsize = 12)
plt.grid(True)
plt.show()
err = np.abs(dfdx - np.cos(x))
print('interior max error %8.2e predicted dx^2/6 %8.2e' % (err[1:-1].max(), dx**2/6))
print('error at x_1 %8.2e predicted dx |f\'\'(x_1)|/2 %8.2e' % (err[0], dx*abs(np.sin(x[0]))/2))
print('error at x_n %8.2e predicted dx |f\'\'(x_n)|/2 %8.2e' % (err[-1], dx*abs(np.sin(x[-1]))/2))
The computed derivative follows $\cos x$ closely, and each printed error is within about $30\%$ of its predicted leading term — the rest is the next term of the expansion, which is not negligible at this coarse spacing. The two end points are about five times worse than the interior here, but the number five is not the point: the interior error is $O(\Delta x^2)$ and the end error is only $O(\Delta x)$, so the gap widens every time the grid is refined. The next cell measures that.
This is the single most common weakness of finite differences on a finite data record, and there are two standard cures. If the data are periodic, wrap around and use the central formula everywhere. Otherwise, use a second-order one-sided formula at the ends, built by the same Taylor bookkeeping as Section 2.4.2 but from the points that do exist:
$$ f'(x_1) \approx \frac{-3f_1 + 4f_2 - f_3}{2\Delta x} + O(\Delta x^2), \qquad f'(x_n) \approx \frac{3f_n - 4f_{n-1} + f_{n-2}}{2\Delta x} + O(\Delta x^2) . $$
# how do the errors shrink as the grid is refined?
print(' n dx interior ends, 1st order ends, 2nd order')
for n in [30, 60, 120, 240]:
x = np.linspace(0.1, 3, n)
dx = x[1] - x[0]
f = np.sin(x)
# first-order ends, as above
d1 = np.zeros(n)
d1[0] = (f[1] - f[0]) / dx # forward diff at f(x_1)
for i in range(1, n - 1):
d1[i] = (f[i + 1] - f[i - 1]) / (2 * dx) # central in between
d1[-1] = (f[-1] - f[-2]) / dx # backward diff at f(x_n)
# second-order ends, same interior
d2 = d1.copy()
d2[0] = (-3*f[0] + 4*f[1] - f[2]) / (2 * dx)
d2[-1] = (3*f[-1] - 4*f[-2] + f[-3]) / (2 * dx)
e1 = np.abs(d1 - np.cos(x))
e2 = np.abs(d2 - np.cos(x))
print('%4d %7.4f %8.2e %8.2e %8.2e'
% (n, dx, e1[1:-1].max(), max(e1[0], e1[-1]), max(e2[0], e2[-1])))
Read the table down each column. Every time $n$ doubles, the interior error falls by a factor of $4$ and so does the second-order end error, while the first-order end error falls by only $2$. By $n = 240$ the cheap end formula is $36$ times worse than the interior, and it alone would decide the accuracy of anything built on top of these derivatives. The second-order ends stay within a factor of two of the interior at every resolution, which is what "matching the interior order" means in practice.
One remark on the code. The interior loop can be replaced by a single vectorized expression, in which the shifted arrays f[2:] and f[:-2] play the roles of $f_{k+1}$ and $f_{k-1}$ — for large grids this matters a great deal. NumPy packages the whole scheme as np.gradient, which uses central differences in the interior and, with edge_order = 2, the same one-sided formulas at the ends.
n = 30
x = np.linspace(0.1, 3, n)
dx = x[1] - x[0]
f = np.sin(x)
# the same thing again, without the loop
dfdx = np.zeros(n)
dfdx[1:-1] = (f[2:] - f[:-2]) / (2 * dx) # central in between
dfdx[0] = (-3*f[0] + 4*f[1] - f[2]) / (2 * dx) # 2nd-order forward
dfdx[-1] = (3*f[-1] - 4*f[-2] + f[-3]) / (2 * dx) # 2nd-order backward
print('max difference from np.gradient %8.1e' % np.max(np.abs(np.gradient(f, dx, edge_order = 2) - dfdx)))
3.1. Why Numerical Integration?¶
3.1.1. Motivation¶
Most integrals that arise in engineering cannot be evaluated by the rules of calculus.
- There is no antiderivative in closed form. Even innocent-looking integrands such as $e^{-x^2}$, $\sin(x)/x$ or $\sqrt{1 + \cos^2 x}$ (the arc length of a sine curve) have none.
- The integrand is sampled data — a force record to be turned into an impulse, a flow rate to be turned into a volume, a power signal to be turned into energy. There is no formula to integrate at all.
- The integrand is the output of a simulation or a black box, available one value at a time.
In each case the only thing we can do is evaluate $f$ at a finite set of points and combine those numbers into an estimate of the area. Every method below is of the form
$$ \int_a^b f(x)\,dx \;\approx\; \sum_k w_k\,f(x_k) , $$
a weighted sum of samples. The methods differ only in where the points $x_k$ are placed and what the weights $w_k$ are.
3.1.2. The Riemann Sum¶
The definition of the integral is already a numerical method waiting to be used. Divide $[a, b]$ into $N$ equal pieces of width
$$ \Delta x = \frac{b - a}{N} , \qquad x_k = a + k\,\Delta x , \qquad k = 0, 1, \dots, N , $$
so that $x_0 = a$ and $x_N = b$. Then
$$ \int_a^b f(x)\,dx = \lim_{N \to \infty} \sum_{k=1}^{N} f\!\left(a + \frac{b-a}{N}k\right)\!\left(\frac{b-a}{N}\right) = \lim_{\Delta x \to 0} \sum_{k=1}^{N} f(x_k)\,\Delta x . $$
Each term $f(x_k)\,\Delta x$ is the area of a rectangle of width $\Delta x$ whose height is the value of $f$ at the right edge of the piece. Taking the height at the left edge instead gives the equally valid
$$ \int_a^b f(x)\,dx = \lim_{\Delta x \to 0} \sum_{k=0}^{N-1} f(x_k)\,\Delta x . $$
As with the derivative, the numerical method is what is left when the limit is dropped and a finite $\Delta x$ is kept. The two choices of edge give the two rectangle rules of Section 3.2, and everything that follows is an attempt to do better than a rectangle.
Throughout we use the test integral
$$ \int_0^{10} \sin x \, dx = \bigl[-\cos x\bigr]_0^{10} = 1 - \cos 10 = 1.8390715\ldots , $$
whose exact value is known, so that every approximation can be compared with the truth.
3.2.1. Left and Right Sums¶
Keeping $\Delta x$ finite in the two limits above gives the left-sided and right-sided rectangle rules,
$$ \int_a^b f(x)\,dx \;\approx\; \sum_{k=0}^{N-1} f(x_k)\,\Delta x \qquad\text{and}\qquad \int_a^b f(x)\,dx \;\approx\; \sum_{k=1}^{N} f(x_k)\,\Delta x . $$
Both are sums of $N$ rectangle areas; they differ only in which endpoint of each piece supplies the height. Note that each uses $N$ of the $N+1$ available samples, the left rule ignoring $f(b)$ and the right rule ignoring $f(a)$.
We take $a = 0$, $b = 10$, $\Delta x = 0.5$, and sample $\sin x$ on that grid. The dashed line in the plot is the staircase whose area the left-sided rule computes.
a = 0
b = 10
dx = 0.5
x = np.arange(a, b + dx, dx)
fx = np.sin(x)
n = len(x) # number of data points
print('number of data points n =', n)
print('number of rectangles N =', n - 1)
exact = 1 - np.cos(10) # the value we are trying to reproduce
print('exact integral ', exact)
# Plot
plt.plot(x, fx, '-k', linewidth = 2)
plt.step(x, fx, '--r', where = 'post') # the staircase of the left-sided rule
plt.grid(True)
plt.xlabel('$x$', fontsize = 12)
plt.show()
The picture below shows the two sets of rectangles themselves. Where $f$ is increasing the left rectangles sit below the curve and the right rectangles above it; where $f$ is decreasing the roles are reversed. That observation is the whole story of the error.
xf = np.linspace(a, b, 400)
plt.figure(figsize = (13, 3.6))
plt.subplot(1, 2, 1)
plt.bar(x[:-1], fx[:-1], width = dx, align = 'edge', color = 'c', edgecolor = 'k', linewidth = 0.5, alpha = 0.5)
plt.plot(xf, np.sin(xf), '-k', linewidth = 2)
# plt.grid(True)
plt.xlim([0, 10])
plt.title('left-sided rectangles')
plt.subplot(1, 2, 2)
plt.bar(x[:-1], fx[1:], width = dx, align = 'edge', color = 'y', edgecolor = 'k', linewidth = 0.5, alpha = 0.5)
plt.plot(xf, np.sin(xf), '-k', linewidth = 2)
plt.xlim([0, 10])
# plt.grid(True)
plt.title('right-sided rectangles')
plt.show()
The two sums are written below exactly as they read: add one rectangle area at a time.
# left-sided rectangle rule
areaL = 0
for i in range(n - 1):
areaL = areaL + fx[i]*dx
# right-sided rectangle rule
areaR = 0
for i in range(n - 1):
areaR = areaR + fx[i + 1]*dx
print('areaL %10.6f error %+9.6f' % (areaL, areaL - exact))
print('areaR %10.6f error %+9.6f' % (areaR, areaR - exact))
A loop is the clearest way to write a sum once, but NumPy can add the whole array at a stroke. The left rule adds the first $N$ samples, the right rule the last $N$, and both multiply by $\Delta x$.
# the same two sums, without the loop
areaL = np.sum(fx[:-1])*dx
areaR = np.sum(fx[1:])*dx
print('areaL %10.6f' % areaL)
print('areaR %10.6f' % areaR)
print('\nareaR - areaL %+9.6f' % (areaR - areaL))
print('( f(b) - f(a) ) * dx %+9.6f' % ((fx[-1] - fx[0])*dx))
3.2.2. Why the Rectangle Rules Are First Order¶
The two sums differ by $0.272$, which is a large disagreement for two approximations of the same number. The last check explains it exactly. The two sums contain the same interior samples $f(x_1), \dots, f(x_{N-1})$; only the two ends differ, so
$$ \text{areaR} - \text{areaL} = \bigl(f(b) - f(a)\bigr)\,\Delta x , $$
which the code confirms to the last digit. The difference shrinks like $\Delta x$, and neither sum is closer to the truth than that difference — so both must carry an error of order $\Delta x$.
The same conclusion follows from a Taylor expansion on one piece. With $x_{k+1} = x_k + \Delta x$,
$$ \int_{x_k}^{x_{k+1}} f(x)\,dx = \Delta x\, f(x_k) + \frac{\Delta x^2}{2}f'(x_k) + O(\Delta x^3) , $$
so a single left rectangle misses by $\tfrac{\Delta x^2}{2}f'(x_k)$. Summing over all $N$ pieces, and recognizing $\Delta x\sum f'(x_k)$ as a Riemann sum for $\int f'$,
$$ \int_a^b f(x)\,dx - \text{areaL} \;\approx\; \frac{\Delta x}{2}\int_a^b f'(x)\,dx = \frac{\Delta x}{2}\bigl(f(b) - f(a)\bigr) . $$
The error of a single rectangle is $O(\Delta x^2)$, but there are $N \propto 1/\Delta x$ of them, so the total error is $O(\Delta x)$ — first-order accurate. Halving $\Delta x$ only halves the error. The right rule has the same error with the opposite sign, which is exactly what the next section exploits.
3.3.1. Averaging the Two Rectangles¶
The left rule misses by $+\tfrac{\Delta x}{2}(f(b) - f(a))$ and the right rule by $-\tfrac{\Delta x}{2}(f(b) - f(a))$. Averaging them must cancel that term. Geometrically, replacing each rectangle by the average of the left and right rectangle is the same as joining $(x_k, f_k)$ to $(x_{k+1}, f_{k+1})$ by a straight line: the shape becomes a trapezoid of area $\tfrac{\Delta x}{2}(f_k + f_{k+1})$, and the rule is
$$ \int_a^b f(x)\,dx \;\approx\; \sum_{k=0}^{N-1} \frac{\Delta x}{2}\bigl(f_k + f_{k+1}\bigr) . $$
Collecting terms, every interior sample appears in two trapezoids and the two end samples in one each, so the rule can also be written with explicit weights,
$$ \int_a^b f(x)\,dx \;\approx\; \Delta x\left(\frac{f_0}{2} + f_1 + f_2 + \cdots + f_{N-1} + \frac{f_N}{2}\right) . $$
The trapezoidal rule therefore costs nothing beyond the rectangle rules: it uses the same $N + 1$ samples and only weights the two ends by one half.
# trapezoidal rule
Z = 0
for i in range(n - 1):
Z = Z + (dx/2)*(fx[i + 1] + fx[i])
print('areaZ %10.6f error %+9.6f' % (Z, Z - exact))
# it really is the average of the two rectangle rules
print('\n(areaL + areaR) / 2 %10.6f' % ((areaL + areaR)/2))
# and NumPy has it built in
print('np.trapezoid(fx, x) %10.6f' % np.trapezoid(fx, x))
The three numbers agree exactly. NumPy's np.trapezoid is the direct counterpart of MATLAB's trapz, and it takes the sampled values and the sample points, so it works on measured data just as well as on a formula. (In NumPy versions before $2.0$ the same function is called np.trapz.)
a = 0
b = 10
dx = 0.5
x = np.arange(a, b + dx, dx)
fx = np.sin(x)
n = len(x) # number of data points
print('number of data points n =', n)
print('number of trapezoids N =', n - 1)
exact = 1 - np.cos(10) # the exact value
print('exact integral ', exact)
# Trapezoidal Rule
integral = np.sum((fx[:-1] + fx[1:]) * dx / 2)
print('trapezoidal integral ', integral)
print('absolute error ', abs(integral - exact))
xf = np.linspace(a, b, 400)
plt.figure(figsize = (13, 3.6))
plt.subplot(1, 2, 1)
# Plot the trapezoids
for i in range(n - 1):
x_trap = [x[i], x[i], x[i+1], x[i+1]]
y_trap = [0, fx[i], fx[i+1], 0]
plt.fill(x_trap, y_trap, color = 'b', alpha = 0.5)
plt.plot([x[i], x[i+1]], [fx[i], fx[i+1]], 'k-', linewidth = 1)
# Exact function
plt.plot(xf, np.sin(xf), '-k', linewidth = 2)
plt.xlim([0, 10])
plt.title('Trapezoidal Rule')
plt.show()
3.3.2. Why the Trapezoidal Rule Is Second Order¶
The error of the trapezoidal rule on the whole interval is
$$ \int_a^b f(x)\,dx - \text{areaZ} \;=\; -\frac{\Delta x^2}{12}\bigl(f'(b) - f'(a)\bigr) + O(\Delta x^4) , $$
so it is second-order accurate: halving $\Delta x$ divides the error by four. The cancellation of the first-order term is the same phenomenon met in Chapter 2, where averaging the forward and backward differences produced the second-order central difference. Symmetry buys an order of accuracy, and it is free.
The formula also predicts the error, not just its order, and at $\Delta x = 0.5$ it does so to three digits.
print('actual error %+9.6f' % (Z - exact))
print('predicted error %+9.6f' % (dx**2/12*(np.cos(10) - np.cos(0))))
Remark (when the trapezoidal rule is extraordinary). The error depends only on the derivatives at the two endpoints. If $f$ is smooth and periodic and $[a, b]$ spans a whole number of periods, then $f'(b) = f'(a)$, every term of the expansion cancels, and the trapezoidal rule converges faster than any power of $\Delta x$. This is why simple equally-spaced sums are the method of choice for periodic problems, and it is the basis of spectral methods.
Every rule so far has placed its sample points on a grid, and every one has become more accurate as that grid was refined. The last section of this chapter abandons the grid altogether.
3.4.1. Sampling Instead of a Grid¶
Monte Carlo integration is a class of algorithms that rely on repeated random sampling to obtain a numerical result. Unlike the rules of Sections 3.2 and 3.3, which evaluate the integrand at fixed, equally spaced points, a Monte Carlo method evaluates it at points drawn at random from the domain of integration.
Named after the casino in Monaco, these methods were developed in the 1940s by Stanislaw Ulam, John von Neumann and Nicholas Metropolis as part of the Manhattan Project. They are now indispensable in physics, financial engineering, statistical mechanics and machine learning.
The formulation begins with the integral of a continuous function $f$ over a domain $\Omega = [a, b]$ whose total length is $|\Omega| = b - a$,
$$ I = \int_{\Omega} f(x)\,dx . $$
If $X$ is a random variable drawn uniformly from $\Omega$, its probability density is the constant $p(x) = 1/|\Omega|$, and the integral can be rewritten as an expected value:
$$ I = |\Omega| \int_{\Omega} f(x)\,\frac{1}{|\Omega|}\,dx = |\Omega| \cdot \mathbb{E}[f(X)] . $$
An expected value is approximated by a sample mean. Drawing $N$ independent points $x_1, \dots, x_N$ uniformly from $\Omega$ gives the Monte Carlo estimate
$$ \int_{\Omega} f(x)\,dx \;\approx\; \frac{1}{N}\sum_{j=1}^{N} f(x_j)\,|\Omega| , \qquad x_j \sim \mathcal{U}(\Omega) . $$
This is still the weighted sum of samples introduced in Section 3.1.1, $\sum_k w_k f(x_k)$. What has changed is that the points are no longer chosen and the weights are all equal to $|\Omega|/N$. Read geometrically, each sample defines a rectangle of width $b - a$ and height $f(x_j)$, and the estimate is the average area of $N$ such rectangles. With $N = 4$,
$$ I \approx \frac{1}{4}\Bigl[ f(x_1)(b-a) + f(x_2)(b-a) + f(x_3)(b-a) + f(x_4)(b-a) \Bigr] , $$
and as $N \to \infty$ the average area of these stochastic rectangles converges to the exact area under the curve.
3.4.2. Hit-or-Miss: the Area of a Circle¶
There is a second, even more elementary way to see the idea, and it applies to a region rather than to a function. Let $\Omega$ be a domain of unknown area $A$, contained in a simple bounding region $B$ whose area $A_B$ is known. If $N$ points are scattered uniformly over $B$, the probability that any one of them lands in $\Omega$ is the ratio of the areas,
$$ P = \frac{\text{Area}(\Omega)}{\text{Area}(B)} = \frac{A}{A_B} . $$
By the law of large numbers, the observed fraction converges to that probability as $N$ grows, so with $N_{\text{inside}}$ points landing in $\Omega$,
$$ \frac{N_{\text{inside}}}{N} \approx \frac{A}{A_B} \qquad\Longrightarrow\qquad A \approx A_B \cdot \frac{N_{\text{inside}}}{N} . $$
This is the hit-or-miss technique. As a demonstration we estimate the area of the unit circle $x^2 + y^2 \le 1$, whose exact value is $\pi$, by enclosing it in the square $[-1, 1] \times [-1, 1]$:
- Define the bounding box. The square has area $A_{\text{box}} = 2 \times 2 = 4$.
- Generate samples. Draw $N$ independent uniform points $(x_i, y_i)$ from $[-1, 1] \times [-1, 1]$.
- Inclusion test. Check whether $x_i^2 + y_i^2 \le 1$.
- Compute the estimate. Take $A \approx 4 \cdot N_{\text{inside}}/N$, and plot a subset of the points to see the sampling distribution inside and outside the boundary.
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
# 1. Setup Parameters
n_samples = 1000000
box_area = 4.0
# 2. Generate Random Points
x = np.random.uniform(-1.0, 1.0, n_samples)
y = np.random.uniform(-1.0, 1.0, n_samples)
# 3. Check Points Inside Circle
distance_squared = x**2 + y**2
inside_mask = distance_squared <= 1.0
n_inside = np.sum(inside_mask)
# 4. Calculate Estimated Area
estimated_area = box_area * (n_inside / n_samples)
exact_area = np.pi
error_percent = abs(estimated_area - exact_area) / exact_area * 100
# 5. Display Results
print(f"Total Samples: {n_samples:,}")
print('Estimated Area: %.5f' % estimated_area)
print('Exact Area (Pi): %.5f' % exact_area)
print('Error: %.3f%%' % error_percent)
# 6. Plot Sample Points (Subset of 5,000 points for clarity)
n_plot = 2000
plt.figure(figsize = (6, 6))
plt.scatter(x[:n_plot][inside_mask[:n_plot]], y[:n_plot][inside_mask[:n_plot]], color = 'crimson', s = 2, label = 'Inside')
plt.scatter(x[:n_plot][~inside_mask[:n_plot]], y[:n_plot][~inside_mask[:n_plot]], color = 'royalblue', s = 2, label = 'Outside')
# Boundary Circle
theta = np.linspace(0, 2 * np.pi, 200)
plt.plot(np.cos(theta), np.sin(theta), 'k--', linewidth = 1.5, label = 'Boundary')
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.gca().set_aspect('equal')
plt.title('Monte Carlo Integration: Unit Circle Area')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc = 'upper right')
plt.grid(True, linestyle = ':', alpha = 0.5)
plt.show()
A million samples locate $\pi$ to three decimal places. That is a modest return for so much arithmetic, and the reason is the subject of Section 3.4.4.
3.4.3. Sample Mean: a One-Dimensional Integral¶
The expected-value formulation of Section 3.4.1 is applied next to a definite integral,
$$ I = \int_{0}^{\pi} \bigl( \sin x + 0.5\cos 3x + 1 \bigr)\,dx . $$
- Analytical solution. $\displaystyle I_{\text{exact}} = \Bigl[-\cos x + \tfrac{1}{6}\sin 3x + x\Bigr]_0^{\pi} = 2 + \pi \approx 5.14159$.
- Domain length. $|\Omega| = b - a = \pi$.
- Computation. The integrand and the area to be computed are plotted first; then $N$ points $x_i \sim \mathcal{U}(0, \pi)$ are drawn and the integral is estimated as the mean height times the domain length.
import numpy as np
import matplotlib.pyplot as plt
# 1. Setup Function and Domain [0, pi]
a, b = 0.0, np.pi
x_plot = np.linspace(a, b, 500)
# f(x) = sin(x) + 0.5*cos(3x) + 1.0
fx_plot = np.sin(x_plot) + 0.5 * np.cos(3 * x_plot) + 1.0
# 2. Plot Function Curve
plt.figure(figsize = (6, 4))
plt.plot(x_plot, fx_plot, 'b-', linewidth = 2, label = r'$f(x) = \sin(x) + 0.5\cos(3x) + 1.0$')
# 3. Fill Integration Area [0, pi]
plt.fill_between(x_plot, fx_plot, color = 'skyblue', alpha = 0.4, label = 'Integration Area')
# 4. Display Settings
plt.axhline(0, color = 'black', linewidth = 0.8, linestyle = '--')
plt.axvline(a, color = 'gray', linestyle = '--')
plt.axvline(b, color = 'gray', linestyle = '--')
plt.title(r'Integrand Function $f(x)$ over $[0, \pi]$')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True, linestyle = ':', alpha = 0.6)
# plt.legend(loc='upper right')
plt.show()
import numpy as np
np.random.seed(0)
# 1. Setup Domain and Parameters
a, b = 0.0, np.pi
domain_length = b - a
n_samples = 1000000
# 2. Generate Uniform Random Samples
x = np.random.uniform(a, b, n_samples)
# 3. Evaluate Integrand: f(x) = sin(x) + 0.5*cos(3x) + 1.0
fx = np.sin(x) + 0.5 * np.cos(3 * x) + 1.0
# 4. Compute Monte Carlo Estimate (Domain Length * Mean Height)
estimated_integral = domain_length * np.mean(fx)
# Analytical exact value: [-cos(x) + (1/6)*sin(3x) + 1.0x] from 0 to pi = 2 + pi approx 5.14159
exact_integral = 2.0 + np.pi
error_percent = abs(estimated_integral - exact_integral) / exact_integral * 100
# 5. Display Results
print(f"Total Samples: {n_samples:,}")
print('Estimated Integral: %.5f' % estimated_integral)
print('Exact Integral: %.5f' % exact_integral)
print('Error: %.3f%%' % error_percent)
There is no grid anywhere in that computation, and no formula for weights: the entire method is "evaluate $f$ at random points and average". Nothing in it refers to the dimension of the domain either, which is the property that matters next.
3.4.4. Convergence and the Curse of Dimensionality¶
The estimate is itself a random variable, so a different set of samples gives a different answer. By the central limit theorem its standard deviation is
$$ \text{Error} = \mathcal{O}\!\left(\frac{1}{\sqrt{N}}\right) , $$
which is a discouraging rate: to gain one more decimal digit, $N$ must be multiplied by one hundred. It is far worse than the $O(\Delta x^2)$ of the trapezoidal rule. Section 3.4.3 showed exactly that: the million random points used there reach a relative error of about $10^{-5}$, which the trapezoidal rule attains on the same integral with about one hundred and sixty subintervals.
The crucial property is not the rate itself but what the rate does not depend on. The exponent $1/2$ is the same in one dimension and in one thousand, because it comes from the statistics of averaging and not from the geometry of a grid.
A grid-based rule has no such immunity. Covering a $d$-dimensional cube with $n$ points per side requires $N = n^d$ evaluations in total, so the spacing is $\Delta x \sim n^{-1} = N^{-1/d}$ and the trapezoidal error $O(\Delta x^2)$ becomes
$$ \mathcal{O}\bigl(N^{-2/d}\bigr) . $$
The comparison is then immediate. In one dimension $N^{-2}$ is vastly better than $N^{-1/2}$, and a grid should be used. The two rates coincide at $d = 4$, and beyond that the grid rule falls behind and keeps falling: at $d = 10$ its error decays as $N^{-1/5}$ while the Monte Carlo error still decays as $N^{-1/2}$. This deterioration of grid methods with dimension is the curse of dimensionality, and escaping it is the reason Monte Carlo integration is used at all.
Remark (where high-dimensional integrals come from). Ten dimensions is not an exotic case. An expected value over a system with ten uncertain parameters, a partition function in statistical mechanics, the price of an option depending on many correlated assets, and the expected loss of a machine learning model over its data distribution are all integrals whose dimension is the number of variables in the problem. For these, sampling is not one option among several; it is the only one.
4.1. From Integration to Differential Equations¶
4.1.1. The Initial Value Problem¶
Almost every dynamical model in engineering can be written as a first-order system
$$ \dot{x} = f(x) , \qquad x(0) = x_0 , $$
where $x$ may be a vector of states and $f$ may be a nonlinear function. A second-order equation is put in this form by naming the velocity as an extra state, as Sections 4.2 and 4.5 will do.
The special case in which $f$ is linear,
$$ \dot{x} = A x , \qquad x(0) = x_0 , $$
is a system of first-order linear differential equations, and it has the closed-form solution $x(t) = e^{At}x_0$. That solution is the exception: for a general nonlinear $f$ no formula exists, and we must compute instead.
What "compute" means here is worth stating precisely. A numerical solver does not return a function. Starting from $x_0$ it iterates, producing
$$ x_0 \;\to\; x_1 \;\to\; x_2 \;\to\; \cdots \;\to\; x_N , \qquad x_k \approx x(t_k), \qquad t_k = k\,\Delta t , $$
a finite list of states sampled along the solution. That list is called a trajectory. Every method in this chapter is a rule for getting from $x_k$ to $x_{k+1}$.
4.1.2. Forward and Backward Euler¶
The simplest rules come straight from the finite differences of two sections ago. Replace $\dot{x}$ by a forward difference and evaluate the right-hand side at the point we already know:
$$ \frac{x_{k+1} - x_k}{\Delta t} \approx \dot{x} = f(x_k) \qquad\Longrightarrow\qquad \;x_{k+1} = x_k + \Delta t\, f(x_k)\; $$
This is forward Euler, also called the explicit Euler method. Everything on the right is known, so the new state is obtained by one evaluation of $f$ and one addition. Geometrically it is the tangent line: follow the slope at the current point for one step.
Replacing $\dot{x}$ by a backward difference instead evaluates the right-hand side at the point we are trying to find:
$$ \frac{x_{k+1} - x_k}{\Delta t} \approx \dot x = f(x_{k+1}) \qquad\Longrightarrow\qquad \;x_{k+1} = x_k + \Delta t\, f(x_{k+1}) $$
This is backward Euler, or the implicit Euler method. The unknown $x_{k+1}$ now appears on both sides, so each step requires solving an equation rather than evaluating a formula. For a general nonlinear $f$ that means a rootfinding problem at every step — which is where Newton's method earns its keep. Section 4.4 explains why anyone would pay that price.
For a linear system the implicit equation can be solved once and for all. With $f(x) = Ax$,
$$ \text{forward:}\quad x_{k+1} = x_k + \Delta t\, A x_k = \bigl(I + \Delta t A\bigr)x_k , $$
$$ \text{backward:}\quad x_{k+1} = x_k + \Delta t\, A x_{k+1} \quad\Longrightarrow\quad x_{k+1} = \bigl(I - \Delta t A\bigr)^{-1} x_k , $$
where $I$ is the identity matrix. Both are just a matrix multiplied by the current state; the difference is which matrix.
4.2.1. The Model¶
A mass $m$ on a spring of stiffness $K$ with a damper of constant $c$ obeys
$$ m\ddot{x} = -Kx - c\dot{x} \qquad\Longrightarrow\qquad \ddot{x} + \frac{c}{m}\dot{x} + \frac{K}{m}x = 0 . $$
Introducing the natural frequency and the damping ratio
$$ \omega_0 = \sqrt{\frac{K}{m}} , \qquad \zeta = \frac{c}{2\sqrt{Km}} , $$
puts it in the standard form of a second-order linear differential equation,
$$ \ddot{x} + 2\zeta\omega_0\dot{x} + \omega_0^2 x = 0 . $$
To solve it numerically we need a first-order system, so we name the velocity $v = \dot{x}$ as a second state:
$$ \begin{aligned} \dot{x} &= v \\ \dot{v} &= -\omega_0^2 x - 2\zeta\omega_0 v \end{aligned} \qquad\Longrightarrow\qquad \frac{d}{dt}\begin{bmatrix} x \\ v \end{bmatrix} = \underbrace{\begin{bmatrix} 0 & 1 \\ -\omega_0^2 & -2\zeta\omega_0 \end{bmatrix}}_{A} \begin{bmatrix} x \\ v \end{bmatrix} . $$
This trick — adding the derivative as a state — turns any $n$-th order equation into a first-order system of $n$ equations, which is why solvers only ever need to handle the form $\dot{x} = f(x)$.
The two numbers $\omega_0$ and $\zeta$ fix the eigenvalues of $A$, and therefore the entire behaviour:
- under-damped, $\zeta < 1$: the system oscillates while decaying, at the damped frequency $\omega_d = \omega_0\sqrt{1 - \zeta^2}$;
- critically damped, $\zeta = 1$: the fastest decay without oscillation;
- over-damped, $\zeta > 1$: two real decay rates, no oscillation.
We take $\omega_0 = 2\pi$ (a one-second period) and $\zeta = 1.75$, an over-damped case, released from rest at $x = 2$.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
w = 2*np.pi
d = 1.75 # damping ratio (zeta)
A = np.array([[0, 1], [-w**2, -2*d*w]])
dt = 0.1 # time step
T = 5 # amount of time to integrate
x0 = np.array([2, 0]) # initial condition
print('A =\n', A)
print('\neigenvalues of A ', np.linalg.eigvals(A))
Both eigenvalues are real and negative, confirming the over-damped case, and the true solution decays to zero. Note how different their sizes are: one mode decays on a time scale of $1/2$ second, the other about ten times faster. Remember that gap; Section 4.4 shows that it is what limits the step size.
4.2.2. Coding the Two Euler Methods¶
Each method is a loop over the same three lines: start from $x_0$, apply the update rule, store the result. For this linear system the update is a single matrix multiplication, so the matrix is formed once outside the loop.
n = int(T/dt) # number of steps
# iterate forward Euler
xF = np.zeros((2, n + 1))
tF = np.zeros(n + 1)
xF[:, 0] = x0
for i in range(n):
tF[i + 1] = (i + 1)*dt
xF[:, i + 1] = (np.eye(2) + A*dt) @ xF[:, i]
# iterate backward Euler
xB = np.zeros((2, n + 1))
tB = np.zeros(n + 1)
xB[:, 0] = x0
for i in range(n):
tB[i + 1] = (i + 1)*dt
xB[:, i + 1] = np.linalg.inv(np.eye(2) - A*dt) @ xB[:, i]
# compute better solution using a built-in solver
sol = solve_ivp(lambda t, y: A @ y, (0, T), x0, t_eval = np.arange(0, T + dt, dt))
t = sol.t
y = sol.y
plt.plot(tB, xB[0, :], 'w') # backward euler
plt.plot(tF, xF[0, :], 'y', label = 'Forward Euler') # forward euler
plt.plot(t, y[0, :], 'k--', label = 'Ground Truth') # built-in solver
# plt.legend(['Forward Euler', 'Backward Euler', 'solve_ivp'], fontsize=12)
plt.legend()
# plt.grid(True)
plt.xlim([0, T])
plt.xlabel('$t$', fontsize = 12)
plt.ylabel('$x$', fontsize = 12)
plt.title('$\\Delta t = %g$' % dt, fontsize = 12)
plt.show()
print('x at t = T')
print(' Forward Euler %10.4f' % xF[0, -1])
print(' Backward Euler %10.4f' % xB[0, -1])
print(' solve_ivp %10.4f' % y[0, -1])
Something is wrong. The true solution has decayed to essentially zero by $t = 5$, and backward Euler and the library solver agree that it has — but forward Euler is still oscillating with an amplitude of about $0.24$, and if the plot is continued the oscillation grows. It is not merely inaccurate; it is qualitatively wrong, producing oscillation in a system that cannot oscillate.
The obvious first response is to shrink the step. It works, and the reason it works is the subject of Section 4.4.
plt.plot(tF, xF[0, :], 'w') # forward euler
dt = 0.05 # half as large
n = int(T/dt)
xF = np.zeros((2, n + 1))
tF = np.arange(n + 1)*dt
xF[:, 0] = x0
for i in range(n):
xF[:, i + 1] = (np.eye(2) + A*dt) @ xF[:, i]
xB = np.zeros((2, n + 1))
xB[:, 0] = x0
for i in range(n):
xB[:, i + 1] = np.linalg.inv(np.eye(2) - A*dt) @ xB[:, i]
sol = solve_ivp(lambda t, y: A @ y, (0, T), x0, t_eval = tF)
# plt.plot(tF, xB[0, :], 'w')
plt.plot(tF, xF[0, :], 'y', label = 'Forward Euler')
plt.plot(sol.t, sol.y[0, :], 'k--', label = 'Ground Truth')
# plt.legend(['Forward Euler', 'Backward Euler', 'solve_ivp'], fontsize=12)
# plt.grid(True)
plt.xlim([0, T])
plt.ylabel('$x$', fontsize = 12)
plt.legend()
plt.xlabel('$t$', fontsize = 12)
plt.title('$\\Delta t = %g$' % dt, fontsize = 12)
plt.show()
print('x at t = T Forward %10.3e Backward %10.3e solve_ivp %10.3e'
% (xF[0, -1], xB[0, -1], sol.y[0, -1]))
With $\Delta t = 0.05$ all three results agree to about $10^{-4}$ and decay to zero. The step size mattered far more than the choice of method.
4.3.1. Local and Global Error¶
How wrong is one Euler step? Write the scalar case $\dot{y} = f(y)$ and expand the exact solution about $t_k$:
$$ y(t_k + \Delta t) = y(t_k) + \Delta t\,\frac{dy}{dt}(t_k) + \frac{\Delta t^2}{2!}\frac{d^2y}{dt^2}(c) , \qquad c \in [t_k, t_{k+1}] , $$
where the last term is the exact Taylor remainder. Since $\dot{y}(t_k) = f(y_k)$, the first two terms are exactly the forward Euler step. Everything the method throws away is the third term, so the local error — the error committed by one step, starting from an exact value — is
$$ \epsilon_{k+1} = y(t_k + \Delta t) - y_{k+1} = \frac{\Delta t^2}{2}\frac{d^2y}{dt^2}(c) = O(\Delta t^2) . $$
That looks encouraging, but it is not the number we care about. Reaching a fixed final time $b$ from $a$ takes $N = (b-a)/\Delta t$ steps, and the local errors accumulate. Adding them up gives the global error
$$ E_N = \sum_{k} \epsilon_k \;\approx\; N \cdot \frac{\Delta t^2}{2}\frac{d^2y}{dt^2} = \frac{b-a}{2}\,\Delta t\,\frac{d^2y}{dt^2} = O(\Delta t) . $$
One power of $\Delta t$ is lost to the number of steps. Forward Euler is therefore a first-order method: halving the step size halves the error at the final time. Backward Euler has the same structure with the opposite sign on its remainder, and is also first order.
Remark (this is the same bookkeeping as the rectangle rule). In Chapter 3 a single rectangle had an $O(\Delta x^2)$ error, there were $N \propto 1/\Delta x$ of them, and the total was $O(\Delta x)$. That is not a coincidence: applying forward Euler to the special case $\dot{y} = f(t)$ is the left-rectangle rule, and backward Euler is the right-rectangle rule. The trapezoidal rule has a differential-equation counterpart too — averaging the two slopes gives $x_{k+1} = x_k + \tfrac{\Delta t}{2}\bigl(f(x_k) + f(x_{k+1})\bigr)$, the second-order Crank-Nicolson method.
4.3.2. Measuring the Order¶
The claim is checked as always: halve the step and watch the error at a fixed final time. We compare against the exact solution $x(T) = e^{AT}x_0$, which for this linear system can be computed to machine precision.
T = 1.0
exact = expm(A*T) @ x0
print('exact x(T) = %.8f\n' % exact[0])
print(' dt FE error ratio BE error ratio')
eF_old = eB_old = 0
for dt in [0.05, 0.025, 0.0125, 0.00625, 0.003125]:
n = int(round(T/dt))
xF = x0.astype(float)
xB = x0.astype(float)
MF = np.eye(2) + A*dt
MB = np.linalg.inv(np.eye(2) - A*dt)
for i in range(n):
xF = MF @ xF
xB = MB @ xB
eF = abs(xF[0] - exact[0])
eB = abs(xB[0] - exact[0])
if eF_old == 0:
print('%9.6f %10.3e - %10.3e -' % (dt, eF, eB))
else:
print('%9.6f %10.3e %6.2f %10.3e %6.2f' % (dt, eF, eF_old/eF, eB, eB_old/eB))
eF_old, eB_old = eF, eB
Both ratios settle on $2$: halving $\Delta t$ halves the error, which is first-order convergence. The two methods are almost equally accurate — their error constants differ by only a few percent. Accuracy is not the reason to prefer one over the other.
4.4.1. The Scalar Test Problem¶
Accuracy asks how far a method drifts from the truth over a fixed interval. Stability asks something cruder and more urgent: does the computed solution stay bounded at all? The failure in Section 4.2.2 was a stability failure, and one scalar problem explains it completely.
Take
$$ \dot{y} = \lambda y , \qquad y(0) = y_0 , $$
whose exact solution $y(t) = e^{\lambda t}y_0$ decays whenever $\operatorname{Re}\lambda < 0$. The two methods give
$$ \text{forward:}\quad y_{k+1} = (1 + \lambda\Delta t)\,y_k , \qquad \text{backward:}\quad y_{k+1} = \frac{1}{1 - \lambda\Delta t}\,y_k . $$
Each step multiplies by a fixed number, so after $N$ steps
$$ y_N = (1 + \lambda\Delta t)^N y_0 , \qquad y_N = \left(\frac{1}{1 - \lambda\Delta t}\right)^{\!N} y_0 . $$
A power $g^N$ stays bounded exactly when $|g| \le 1$. The two methods are therefore unstable when
$$ \bigl|1 + \lambda\Delta t\bigr| > 1 \qquad\text{and}\qquad \left|\frac{1}{1 - \lambda\Delta t}\right| > 1 . $$
Write $z = \lambda\Delta t$, a complex number combining the physics ($\lambda$) with our choice ($\Delta t$). Forward Euler is stable for $|1 + z| \le 1$: a disc of radius one centred at $z = -1$. Backward Euler is stable for $|1 - z| \ge 1$: everything outside the disc of radius one centred at $z = +1$ — which includes the entire left half-plane.
That is the whole story. If the true solution decays, backward Euler will decay too, whatever step size is chosen. Forward Euler will only decay if $\lambda\Delta t$ happens to land in its small disc. For a real negative $\lambda$ the condition reduces to
$$ \Delta t \le \frac{2}{|\lambda|} , $$
a hard ceiling set by the fastest decaying mode in the problem, however uninteresting that mode may be.
# the two stability regions in the complex plane z = lambda*dt
re = np.linspace(-3, 3, 400)
im = np.linspace(-3, 3, 400)
X, Y = np.meshgrid(re, im)
Z = X + 1j*Y
plt.figure(figsize = (11, 4.5))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, np.abs(1 + Z), levels = [0, 1], colors = ['c'])
plt.axhline(0, color = 'k', linewidth = 0.5)
plt.axvline(0, color = 'k', linewidth = 0.5)
plt.grid(True)
plt.gca().set_aspect('equal', adjustable = 'box')
plt.xlabel('Re $z$', fontsize = 12)
plt.ylabel('Im $z$', fontsize = 12)
plt.title('Forward Euler: stable where $|1 + z| \\leq 1$', fontsize = 11)
plt.subplot(1, 2, 2)
plt.contourf(X, Y, np.abs(1/(1 - Z)), levels = [0, 1], colors = ['c'])
plt.axhline(0, color = 'k', linewidth = 0.5)
plt.axvline(0, color = 'k', linewidth = 0.5)
plt.grid(True)
plt.gca().set_aspect('equal', adjustable = 'box')
plt.xlabel('Re $z$', fontsize = 12)
plt.title('Backward Euler: stable where $|1 - z| \\geq 1$', fontsize = 11)
plt.show()
The shaded regions are where each method is stable. Forward Euler owns one small disc; backward Euler owns everything except one small disc, and in particular the whole left half-plane where all decaying physics lives. This is what "better stability" means, and it is the reason for tolerating the cost of solving an equation at every step.
4.4.2. Systems: the Eigenvalue Test¶
For a system $\dot{x} = Ax$ the same argument applies to each mode. After $N$ steps,
$$ x_N = \bigl(I + \Delta t A\bigr)^{N} x_0 , \qquad x_N = \bigl(I - \Delta t A\bigr)^{-N} x_0 , $$
and a matrix power stays bounded exactly when every eigenvalue of that matrix has magnitude at most one. So
$$ \text{forward Euler is stable when } \bigl|\text{eig}\bigl(I + \Delta t A\bigr)\bigr| \le 1 , \qquad \text{backward Euler when } \bigl|\text{eig}\bigl((I - \Delta t A)^{-1}\bigr)\bigr| \le 1 . $$
Now we can diagnose Section 4.2.2 exactly. Recall the two eigenvalues of $A$, $\lambda \approx -1.97$ and $\lambda \approx -20.0$. The fast one demands $\Delta t \le 2/20.0 = 0.0999$ — and we used $\Delta t = 0.1$, a hair over the line.
print('eigenvalues of A ', np.linalg.eigvals(A))
print('forward Euler needs dt <= 2/|lambda|_max = %.5f\n' % (2/max(abs(np.linalg.eigvals(A)))))
print(' dt max |eig(I + dt A)| max |eig((I - dt A)^-1)|')
for dt in [0.2, 0.1, 0.05, 0.01]:
gF = max(abs(np.linalg.eigvals(np.eye(2) + A*dt)))
gB = max(abs(np.linalg.eigvals(np.linalg.inv(np.eye(2) - A*dt))))
print('%6.2f %8.4f %s %8.4f %s'
% (dt, gF, 'UNSTABLE' if gF > 1 else 'stable ',
gB, 'UNSTABLE' if gB > 1 else 'stable'))
At $\Delta t = 0.1$ the forward-Euler amplification factor is $1.0019$ — barely above one, which is why the solution in Section 4.2.2 grew slowly rather than exploding. At $\Delta t = 0.2$ it is far above one. Backward Euler stays below one at every step size tested, as the stability region promised.
The same test can be run over a range of damping ratios, and it is worth reading the result carefully.
print(' zeta eig(A) max |eig(I + dt A)| at dt = 0.1')
for dd in [0.25, 1.0, 1.75, 2.0, 20.0]:
Ad = np.array([[0, 1], [-w**2, -2*dd*w]])
ev = np.linalg.eigvals(Ad)
gF = max(abs(np.linalg.eigvals(np.eye(2) + Ad*0.1)))
print('%7.2f %-28s %10.4f %s'
% (dd, np.array2string(ev, precision = 2), gF, 'UNSTABLE' if gF > 1 else 'stable'))
Only one row is stable, and it is not the one a physicist would guess. There are two distinct ways to fall outside the small disc. The lightly damped case $\zeta = 0.25$ has complex eigenvalues lying close to the imaginary axis — and the forward-Euler disc touches that axis only at the origin, so an oscillator with little damping is almost always outside it. At the other end, raising $\zeta$ above one makes the slow mode slower and the fast mode very much faster, until at $\zeta = 20$ the eigenvalues differ by a factor of sixteen hundred and the amplification factor is $24$.
A system whose eigenvalues differ by orders of magnitude like this is called stiff, and stiffness is the single most common reason for an explicit solver to crawl: the step size is dictated by the fastest mode even after that mode has decayed away and contributes nothing to the answer. A physically tamer system can therefore be numerically far harder. Implicit methods exist precisely for these problems, and general-purpose libraries offer both families — solve_ivp defaults to the explicit RK45, and switches to implicit methods such as Radau or BDF when asked.
4.5.1. The Equation of Motion¶
Nothing so far required $f$ to be linear; the matrix form was a convenience. The pendulum shows the general case.
A point mass $m$ on a rigid massless rod of length $L$, with $\theta$ measured from the downward vertical, has kinetic and potential energy
$$ T = \tfrac12 m L^2 \dot\theta^2 , \qquad V = m g L\,(1 - \cos\theta) , $$
and the Euler-Lagrange equation applied to $\mathcal{L} = T - V$ gives
$$ \frac{d}{dt}\frac{\partial \mathcal{L}}{\partial\dot\theta} - \frac{\partial \mathcal{L}}{\partial\theta} = 0 \qquad\Longrightarrow\qquad \ddot\theta = \frac{g}{L}\sin\theta , $$
in which $g$ carries the sign of the gravitational acceleration along the coordinate direction — it is entered as a negative number, so that the equation reads $\ddot\theta = -|g|/L \cdot \sin\theta$ and the hanging position $\theta = 0$ is the stable one. Adding linear damping $-d\,\dot\theta$ and introducing $v = \dot\theta$ gives the first-order system
$$ \begin{aligned} \dot\theta &= v \\ \dot v &= \frac{g}{L}\sin\theta - d\,v \end{aligned} \qquad\Longrightarrow\qquad \frac{d}{dt}\begin{bmatrix}\theta \\ v\end{bmatrix} = f\!\left(\begin{bmatrix}\theta \\ v\end{bmatrix}\right) . $$
The $\sin\theta$ makes this nonlinear, so there is no matrix $A$ and no closed-form solution. For small angles $\sin\theta \approx \theta$ and the equation reduces to the spring-mass-damper of Section 4.2 — which is exactly why that example is worth understanding first.
In code, $f$ becomes a function of $t$ and the state vector $y = [\theta, v]$, returning the list of derivatives. This is the interface every solver expects.
def pendulum_ode(t, y, g, L, d):
return [y[1],
(g/L)*np.sin(y[0]) - d*y[1]]
4.5.2. Simulating It¶
With $g = -10$ and $L = 10$ the ratio $g/L = -1$, so the small-angle period is $2\pi\sqrt{L/|g|} = 2\pi$ seconds. The pendulum is released from $\theta = \pi/4$ at rest and followed for $50$ seconds, with light damping $d = 0.1$.
# define physical constants
g = -10
L = 10
d = 0.1
# define time slices and initial condition
t = np.linspace(0, 50, 500)
y0 = [np.pi/4, 0]
pendulum_sol = solve_ivp(pendulum_ode, (0, 50), y0, t_eval = t, args = (g, L, d,))
y = pendulum_sol.y
print('theta(0) = %7.4f' % y[0, 0])
print('theta(50) = %7.4f' % y[0, -1])
print('small-angle period 2*pi*sqrt(L/|g|) = %.4f s' % (2*np.pi*np.sqrt(L/abs(g))))
plt.figure(figsize = (12, 4))
plt.subplot(1, 2, 1)
plt.plot(t, y[0, :], 'k', linewidth = 1.2)
plt.grid(True)
plt.xlabel('$t$', fontsize = 12)
plt.ylabel('$\\theta$', fontsize = 12)
plt.title('angle against time', fontsize = 11)
plt.subplot(1, 2, 2)
plt.plot(y[0, :], y[1, :], 'k', linewidth = 1.2)
plt.grid(True)
plt.xlabel('$\\theta$', fontsize = 12)
plt.ylabel('$v$', fontsize = 12)
plt.title('phase portrait', fontsize = 11)
plt.show()
# make 3D plot
fig = plt.figure(figsize = (6, 4.5))
ax = fig.add_subplot(111, projection = '3d')
ax.plot(t, y[0, :], y[1, :], 'k', linewidth = 1.0)
ax.set_xlabel('$t$')
ax.set_ylabel('$\\theta$')
ax.set_zlabel('$v$')
plt.show()
The left panel shows the familiar decaying oscillation. The right panel is the phase portrait, the trajectory drawn in the $(\theta, v)$ plane with time eliminated: it spirals inward to the origin, which is the hanging equilibrium. The three-dimensional plot shows how the two pictures are related — the phase portrait is the shadow of the trajectory looking down the time axis.
4.5.3. Forward Euler on a Nonlinear Problem¶
The Euler loop needs no modification for a nonlinear problem. There is no matrix to form, so $f$ is simply called at each step.
dt = 0.1
n = int(50/dt)
yF = np.zeros((2, n + 1))
tF = np.arange(n + 1)*dt
yF[:, 0] = y0
for i in range(n):
f = np.array(pendulum_ode(tF[i], yF[:, i], g, L, d))
yF[:, i + 1] = yF[:, i] + dt*f
plt.plot(tF, yF[0, :], 'k', linewidth = 1.2)
plt.plot(t, y[0, :], 'r--', linewidth = 1.2)
plt.legend(['Forward Euler, $\\Delta t = 0.1$', 'solve_ivp'], fontsize = 11)
plt.grid(True)
plt.xlabel('$t$', fontsize = 12)
plt.ylabel('$\\theta$', fontsize = 12)
plt.show()
print('theta(50) Forward Euler %8.4f solve_ivp %8.4f\n' % (yF[0, -1], y[0, -1]))
# amplitude over the last ten seconds
print('amplitude for t in [40, 50]')
print(' Forward Euler %.4f' % np.abs(yF[0, tF >= 40]).max())
print(' solve_ivp %.4f' % np.abs(y[0, t >= 40]).max())
# linearize about theta = 0 and apply the eigenvalue test
Alin = np.array([[0, 1], [g/L, -d]])
print('\nlinearized eigenvalues ', np.linalg.eigvals(Alin))
print('max |eig(I + dt Alin)| %.10f' % max(abs(np.linalg.eigvals(np.eye(2) + Alin*dt))))
The two curves agree at first and then part company completely. Over the last ten seconds the reference solution has decayed to an amplitude of about $0.10$, while forward Euler is still swinging at $0.68$ — seven times too large. The method has not merely accumulated first-order error; it has lost the damping altogether.
The eigenvalue test says why. Linearizing about $\theta = 0$ gives $\lambda = -d/2 \pm i\sqrt{\omega_0^2 - d^2/4}$ with $\omega_0^2 = |g|/L$, and a line of algebra gives the forward-Euler amplification factor
$$ \bigl|1 + \lambda\Delta t\bigr|^2 = 1 - d\,\Delta t + \omega_0^2\,\Delta t^2 . $$
The middle term is the physical damping, the last term is the energy the method injects, and with $d = 0.1$, $\Delta t = 0.1$ and $\omega_0 = 1$ they cancel exactly — the printed factor is $1$ to machine precision. Forward Euler is reporting a pendulum that never stops swinging, because its own numerical energy production happens to match the damper.
Remove the damping and only the injection is left. An undamped pendulum conserves energy and its linearized eigenvalues are purely imaginary, $\lambda = \pm i\omega_0$; the formula above becomes $|1 + \lambda\Delta t|^2 = 1 + \omega_0^2\Delta t^2 > 1$ for every step size. Forward Euler cannot help adding energy to an oscillator.
dt = 0.05
n = int(50/dt)
yE = np.zeros((2, n + 1))
tE = np.arange(n + 1)*dt
yE[:, 0] = y0
for i in range(n):
f = np.array(pendulum_ode(tE[i], yE[:, i], g, L, 0)) # d = 0, no damping
yE[:, i + 1] = yE[:, i] + dt*f
energy = 0.5*yE[1, :]**2 + (abs(g)/L)*(1 - np.cos(yE[0, :]))
plt.figure(figsize = (12, 4))
plt.subplot(1, 2, 1)
plt.plot(tE, yE[0, :], 'k', linewidth = 1.0)
plt.grid(True)
plt.xlabel('$t$', fontsize = 12)
plt.ylabel('$\\theta$', fontsize = 12)
plt.title('undamped pendulum, Forward Euler', fontsize = 11)
plt.subplot(1, 2, 2)
plt.plot(tE, energy, 'k', linewidth = 1.2)
plt.grid(True)
plt.xlabel('$t$', fontsize = 12)
plt.ylabel('energy per unit mass', fontsize = 12)
plt.title('energy should be constant', fontsize = 11)
plt.show()
print('energy at t = 0 %.5f' % energy[0])
print('energy at t = 50 %.5f (a factor of %.1f)' % (energy[-1], energy[-1]/energy[0]))
The amplitude grows steadily and the energy — which the physics requires to be constant — rises by a factor of nearly six over fifty seconds. No step size removes this; it only slows it down. A problem with this structure calls for a method whose stability region contains part of the imaginary axis, such as a Runge-Kutta method, or a symplectic integrator designed to conserve energy exactly.
5.1. Two-Point Boundary Value Problems¶
5.1.1. Initial Value or Boundary Value¶
Chapter 4 solved the initial value problem: the entire state is known at $t = 0$, and the differential equation says how to march forward. Every method there was a rule for getting from $x_k$ to $x_{k+1}$, and the solution was built one step at a time.
Many problems in engineering are not posed that way. Consider the steady temperature in a rod whose two ends are held at fixed temperatures, or the deflection of a beam that is pinned at both ends. The governing equation is still a second-order differential equation, but the two conditions needed to pin down its solution are given at opposite ends of the interval:
$$ u''(x) = f(x), \qquad a < x < b, \qquad u(a) = \alpha, \quad u(b) = \beta . $$
This is a two-point boundary value problem. Nothing can be marched: at $x = a$ we know $u$ but not $u'$, so there is no complete state to step forward from. The value of the solution in the middle of the interval depends on both ends at once.
That single difference changes the numerical strategy completely.
- In an initial value problem the unknowns are computed one after another, each from the ones before it.
- In a boundary value problem the unknowns are all coupled, and must be computed simultaneously.
A boundary value problem therefore does not turn into a loop; it turns into a system of algebraic equations. And the way to turn a differential equation into algebraic equations is already in hand: replace the derivatives by the finite differences of Chapter 2.
5.1.2. The Model Problem¶
We take the simplest representative problem,
$$ u''(x) = f(x), \qquad 0 < x < 1, \qquad u(0) = \alpha, \quad u(1) = \beta , $$
which describes the steady-state temperature $u$ along a rod of unit length and unit conductivity, heated internally at the rate $-f(x)$, with the two ends held at the prescribed temperatures $\alpha$ and $\beta$. The same equation describes the small deflection of a loaded string and the electrostatic potential of a given charge distribution, and the numerical treatment is identical in all three cases.
For a test we choose a problem whose answer is known,
$$ f(x) = -\pi^2 \sin(\pi x), \qquad \alpha = 0, \quad \beta = 1 , \qquad\Longrightarrow\qquad u(x) = \sin(\pi x) + x , $$
as is verified by differentiating twice and checking the two end values. Everything below is written for a general $f$, $\alpha$ and $\beta$; the exact solution is used only to measure the error.
# the rod: u'' = f(x), u(0) = alpha, u(1) = beta
def f(x):
return -np.pi**2*np.sin(np.pi*x)
def uexact(x): # the exact solution, used only to check the error
return np.sin(np.pi*x) + x
alpha = 0
beta = 1
5.2.1. From Derivatives to a Matrix¶
Lay down the same uniform grid as in Section 2.5,
$$ x_k = k\,\Delta x, \qquad k = 0, 1, \dots, n , \qquad \Delta x = \frac{1}{n} , $$
and write $u_k$ for the approximation to $u(x_k)$. The two end values are not unknown at all: the boundary conditions give $u_0 = \alpha$ and $u_n = \beta$ directly.
$$ \begin{aligned} k = 0: \quad & u_0 = \alpha \\ k = n: \quad & u_{n} = \beta \end{aligned} $$
The unknowns are the $n - 1$ interior values $u_1, \dots, u_{n-1}$. At each interior point, replace $u''$ by the three-point central formula of Section 2.4.1:
$$ \frac{u_{k-1} - 2u_k + u_{k+1}}{\Delta x^2} = f(x_k), \qquad k = 1, 2, \dots, n-1 . $$
That is one equation per interior point, and there are exactly as many equations as unknowns. The differential equation has become a linear system. To see its structure, write the equations out one by one, multiplied through by $\Delta x^2$ for the moment so that the pattern is easy to read:
$$ \begin{aligned} k = 1: \quad & u_0 - 2u_1 + u_2 = \Delta x^2 f(x_1), \\ k = 2: \quad & u_1 - 2u_2 + u_3 = \Delta x^2 f(x_2), \\ & \;\;\vdots \\ k = n-2: \quad & u_{n-3} - 2u_{n-2} + u_{n-1} = \Delta x^2 f(x_{n-2}), \\ k = n-1: \quad & u_{n-2} - 2u_{n-1} + u_n = \Delta x^2 f(x_{n-1}). \end{aligned} $$
Each equation involves exactly three consecutive grid values, the point itself and its two neighbours, with the coefficients $1, -2, 1$.
Two of these equations are special. The one at $k = 1$ contains $u_0$, and the one at $k = n-1$ contains $u_n$; both of those are known numbers, $\alpha$ and $\beta$. A known quantity does not belong on the side of the unknowns, so it is subtracted from both sides, and the first and last equations become
$$ \begin{aligned} k = 1: \quad & -2u_1 + u_2 = \Delta x^2 f(x_1) - \alpha, \\ k = n-1: \quad & u_{n-2} - 2u_{n-1} = \Delta x^2 f(x_{n-1}) - \beta . \end{aligned} $$
The interior equations, $k = 2, \dots, n-2$, are unchanged. Now every equation contains only the unknowns $u_1, \dots, u_{n-1}$ on the left, and only known quantities on the right.
The system is assembled into a matrix in the usual way: the $k$-th equation becomes the $k$-th row, and the coefficient of $u_j$ in that equation goes into column $j$. From the list above, row $k$ therefore has the entry $1$ in column $k-1$, the entry $-2$ in column $k$, and the entry $1$ in column $k+1$, and all other entries are zero. In the first row the entry in column $0$ is missing, and in the last row the entry in column $n$ is missing, because the values they would multiply, $u_0$ and $u_n$, have been moved to the right-hand side. Dividing by $\Delta x^2$ again and writing the system as $Au = b$ with $u = (u_1, \dots, u_{n-1})^T$,
$$ \underbrace{ \frac{1}{\Delta x^2} \begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ & & & 1 & -2 \end{bmatrix} }_{A} \underbrace{ \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_{n-2} \\ u_{n-1} \end{bmatrix} }_{u} = \underbrace{ \begin{bmatrix} f(x_1) - \alpha/\Delta x^2 \\ f(x_2) \\ \vdots \\ f(x_{n-2}) \\ f(x_{n-1}) - \beta/\Delta x^2 \end{bmatrix} }_{b} . $$
Reading the matrix row by row against the equations above shows where every entry comes from: the diagonal $-2$ is the coefficient of $u_k$ in the $k$-th equation, the two $1$'s beside it are the coefficients of $u_{k-1}$ and $u_{k+1}$, and the $\alpha$ and $\beta$ in the right-hand side are the boundary values that left the first and last equations.
The matrix is the discrete Laplacian $(1, -2, 1)/\Delta x^2$ met in Section 2.4.1, now written out for a whole grid. It is tridiagonal: every row contains at most three nonzero entries, because the stencil couples each point only to its two immediate neighbours. The boundary values appear nowhere in the matrix and only in the first and last entries of the right-hand side: changing $\alpha$ or $\beta$ changes $b$, not $A$.
This is the difference from Chapter 4 in one sentence: a time-stepping method applies the stencil sequentially, computing each new value from values already known, whereas a boundary value method assembles it into a matrix and solves once, because no interior value can be computed before the others.
n = 20 # number of subintervals
x = np.linspace(0, 1, n + 1)
dx = x[1] - x[0]
print('grid spacing dx = %.4f' % dx)
print('unknowns n - 1 = %d' % (n - 1))
5.2.2. Building and Solving the System¶
The matrix is built one row at a time, which is the literal transcription of "one finite difference equation per interior point". The two if statements are there only because the first and last rows have no left or no right neighbour.
A = np.zeros((n - 1, n - 1))
# one row of A is one finite difference equation
for i in range(n - 1):
A[i, i] = -2
if i > 0:
A[i, i - 1] = 1
if i < n - 2:
A[i, i + 1] = 1
A = A/dx**2
print('top left corner of A, in units of 1/dx^2')
print(A[0:4, 0:4]*dx**2)
The right-hand side is the source term sampled at the interior points, with the two boundary values subtracted at the two ends. One call to a linear solver then produces all $n - 1$ interior values at once, and the two known end values are put back in place.
b = np.zeros(n - 1)
for i in range(n - 1):
b[i] = f(x[i + 1])
# the two known boundary values move to the right-hand side
b[0] = b[0] - alpha/dx**2
b[-1] = b[-1] - beta/dx**2
u = np.zeros(n + 1)
u[0] = alpha
u[-1] = beta
u[1:-1] = np.linalg.solve(A, b)
print('max error %9.2e' % np.max(np.abs(u - uexact(x))))
# Plot
xf = np.linspace(0, 1, 400)
plt.plot(xf, uexact(xf), 'k--')
plt.plot(x, u, 'ro', markersize = 5)
plt.xlim(0, 1)
plt.ylim(0, 2)
plt.xlabel('x', fontsize = 12)
plt.ylabel('u(x)', fontsize = 12)
plt.legend(['exact solution', 'finite difference'], fontsize = 12)
plt.show()
Twenty subintervals already reproduce the solution to about two parts in a thousand, and the computed points are indistinguishable from the exact curve by eye. Note what was not needed: no initial slope, no time stepping, no stability restriction on the step size. The boundary value problem is not integrated, it is solved.
This chapter closes the circle opened in Chapter 2. A derivative was replaced by a difference of neighbouring values; that difference, applied at every point of a grid and collected into a matrix, converts a differential equation into a linear system. The finite difference method for partial differential equations is this same construction carried out in more than one variable.
6.1. Introduction to Partial Differential Equations¶
6.1.1. Two Independent Variables¶
Every problem so far has had one independent variable: time in Chapter 4, position in Chapter 5. A partial differential equation involves a function of several variables, here $u(x, t)$ of space and time, and its partial derivatives. The three classical examples are
$$ \begin{align*} u_t &= \alpha\,u_{xx} \quad \text{(heat equation, diffusion)}, \\\\ u_{tt} &= c^2\,u_{xx} \quad \text{(wave equation)}, \\\\ u_{xx} + u_{yy} &= 0 \quad \text{(Laplace equation, equilibrium)} \end{align*} $$
Chapter 5 was already the last kind in one dimension: with no time in the problem, $u'' = f$ described equilibrium. This chapter takes the heat equation and puts the time back in.
Remark. We discretize with finite differences, replacing each derivative by a difference quotient, as in Chapters 2 and 5. The finite element and finite volume methods, standard in structural mechanics and fluid dynamics, reach a similar system of algebraic equations by a different route.
6.1.2. The One-Dimensional Heat Equation¶
The model problem is the transient temperature of a thin rod of length $L$,
$$ \frac{\partial u}{\partial t} = \alpha\,\frac{\partial^2 u}{\partial x^2} , \qquad 0 < x < L , \qquad t > 0 , $$
where $\alpha > 0$ is the thermal diffusivity. The equation says that a point heats up in proportion to how much colder it is than the average of its two neighbours — which is exactly what the second difference of Section 2.4.1 measures.
Because both a time derivative and a space derivative appear, the problem needs both kinds of data met so far: an initial condition $u(x, 0) = u_0(x)$ as in Chapter 4, and boundary conditions at the two ends for all $t$ as in Chapter 5. Here both ends are held at zero, $u(0, t) = u(L, t) = 0$ — an ice bath at each end.
Discretize both variables,
$$ x_i = i\,\Delta x \quad (i = 0, \dots, N_x) , \qquad t_n = n\,\Delta t \quad (n = 0, 1, 2, \dots) , \qquad u_i^n \approx u(x_i, t_n) , $$
the subscript for space and the superscript for time. The computed solution is no longer a curve but a table: one row per time level, each row a snapshot of the whole rod.
For a test we take the initial temperature $u_0(x) = \sin(\pi x / L)$, which is an eigenfunction of $\partial^2/\partial x^2$ and therefore decays in place without changing shape:
$$ u(x, t) = e^{-\alpha \pi^2 t / L^2}\,\sin\!\left(\frac{\pi x}{L}\right) . $$
Every approximation below can be compared against this.
# the rod: u_t = alpha u_xx, u(0,t) = u(L,t) = 0
L = 1.0 # length of the rod
alpha = 1.0 # thermal diffusivity
T = 0.1 # total simulation time
def uexact(x, t): # used only to check the error
return np.exp(-alpha*np.pi**2*t/L**2)*np.sin(np.pi*x/L)
Nx = 50 # number of spatial segments
dx = L/Nx
x = np.linspace(0, L, Nx + 1)
6.2.1. Forward in Time, Central in Space¶
Replace each derivative by the formula already derived for it. For the time derivative use the forward difference of Section 2.2, which is the forward Euler method of Chapter 4; for the space derivative use the three-point central stencil of Section 2.4.1, the same one that built the matrix in Chapter 5:
$$ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha\left(\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\right) . $$
Only one quantity refers to the new time level, so it can be solved for directly:
$$ u_i^{n+1} = u_i^n + r\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right) , \qquad r = \frac{\alpha\,\Delta t}{\Delta x^2} . $$
This is the explicit method, or FTCS (forward time, central space). Each new value is a weighted average of three old ones, no system has to be solved, and the whole method is one line inside a loop. The single dimensionless mesh ratio $r$ contains everything about the discretization, and the entire behaviour of the scheme turns out to depend on it alone.
The update formula is one line, but it is applied once for every interior point. Writing it out for $i = 1, 2, 3, \dots$ at a fixed time level $n$ shows what a single sweep of the loop does:
$$ \begin{aligned} i = 1: \quad u_1^{n+1} &= u_1^n + r\left(u_2^n - 2u_1^n + u_0^n\right), \\ i = 2: \quad u_2^{n+1} &= u_2^n + r\left(u_3^n - 2u_2^n + u_1^n\right), \\ i = 3: \quad u_3^{n+1} &= u_3^n + r\left(u_4^n - 2u_3^n + u_2^n\right), \\ &\;\;\vdots \\ i = N_x - 1: \quad u_{N_x-1}^{n+1} &= u_{N_x-1}^n + r\left(u_{N_x}^n - 2u_{N_x-1}^n + u_{N_x-2}^n\right). \end{aligned} $$
Three things are visible in this list that the compact formula hides.
Every line has the same shape, and only the subscript moves by one from line to line. That is the inner loop, for i in range(1, N_x), with the subscript as the array index.
Every quantity on the right-hand side carries the superscript $n$, and only the quantity on the left carries $n+1$. The whole right-hand side is therefore known before the sweep begins, and the lines can be evaluated in any order; none of them needs a result from another line. This is why the method is called explicit, and why the code must read from the old array and write into a new one rather than overwrite in place.
The first line uses $u_0^n$ and the last uses $u_{N_x}^n$. These are the boundary values, which are never computed by the formula but supplied by the boundary conditions; the loop runs over $i = 1, \dots, N_x - 1$ and leaves the two ends alone.
Once the sweep is complete, the new level replaces the old one, $n \leftarrow n + 1$, and the same list of lines is evaluated again. That is the outer loop.
r = 0.45 # mesh ratio
dt = r*dx**2/alpha
Nt = int(round(T/dt))
print('dx = %.4f dt = %.2e' % (dx, dt))
print('r = %.2f number of time steps = %d' % (r, Nt))
print(Nt)
print(Nx)
print(range(Nt))
print(range(1, Nx))
u = np.sin(np.pi*x/L) # initial temperature
u[0] = 0 # both ends held at zero
u[-1] = 0
unew = np.zeros(Nx + 1)
plt.plot(x, u, 'k', linewidth = 2)
t = 0
for n in range(Nt):
# one explicit update at every interior point
for i in range(1, Nx):
unew[i] = u[i] + r*(u[i+1] - 2*u[i] + u[i-1])
unew[0] = 0 # boundary conditions
unew[-1] = 0
u = unew.copy()
t = t + dt
for tp in [0.01, 0.03, 0.05, 0.10]:
if abs(t - tp) < dt/2:
plt.plot(x, u, linewidth = 1.5)
plt.xlim([0, 1])
plt.xlabel('x', fontsize = 12)
plt.ylabel('u(x, t)', fontsize = 12)
plt.legend(['t = 0', 't = 0.01', 't = 0.03', 't = 0.05', 't = 0.10'], fontsize = 9, loc = 'upper right')
plt.show()
print('max error at t = %.2f %8.2e' % (t, np.max(np.abs(u - uexact(x, t)))))
The profile keeps its shape and decays, as the single exponential predicts, and after five hundred and fifty-six steps it is still correct to two parts in ten thousand. The method is doing what it should.
# the whole history in one array: one row per time level, one column per grid point
U = np.zeros((Nt + 1, Nx + 1))
U[0, :] = np.sin(np.pi*x/L) # initial temperature
for n in range(Nt):
for i in range(1, Nx):
U[n+1, i] = U[n, i] + r*(U[n, i+1] - 2*U[n, i] + U[n, i-1])
U[n+1, 0] = 0 # boundary conditions
U[n+1, -1] = 0
tgrid = np.arange(Nt + 1)*dt
X, Tg = np.meshgrid(x, tgrid) # X[n, i] = x_i, Tg[n, i] = t_n, U[n, i] = u(x_i, t_n)
plt.figure(figsize = (9, 6))
ax = plt.axes(projection = '3d')
ax.plot_surface(X, Tg, U, cmap = 'viridis', rstride = 10, cstride = 1, linewidth = 0, alpha = 0.5)
for tp in [0, 0.01, 0.03, 0.05, 0.10]:
n = int(round(tp/dt)) # time level closest to tp
ax.plot(x, tp*np.ones(Nx + 1), U[n, :], linewidth = 2, color = 'r')
ax.set_xlabel('x', fontsize = 12)
ax.set_ylabel('t', fontsize = 12)
ax.set_zlabel('u(x, t)', fontsize = 12)
ax.view_init(elev = 30, azim = -45)
plt.show()
6.2.2. The Stability Condition¶
The choice $r = 0.45$ was not arbitrary. Group the update rule differently,
$$ u_i^{n+1} = (1 - 2r)\,u_i^n + r\left(u_{i+1}^n + u_{i-1}^n\right) , $$
and read it as a weighted average of the three old values, with weights $r$, $1 - 2r$, $r$ that add up to one. As long as $r \le 1/2$ all three weights are non-negative, so the new value lies between the smallest and the largest of its neighbours and nothing can grow. As soon as $r > 1/2$ the centre weight turns negative: a point that is hotter than its neighbours is pushed below them, overshooting further at every step. The result is the stability condition
$$ r = \frac{\alpha\,\Delta t}{\Delta x^2} \le \frac{1}{2} \qquad\Longleftrightarrow\qquad \Delta t \le \frac{\Delta x^2}{2\alpha} . $$
The same conclusion follows from the eigenvalue test of Section 4.4.2. In matrix form one explicit step is $u^{n+1} = (I + rL)u^n$ with $L$ the discrete Laplacian of Chapter 5, whose eigenvectors are the grid sines and whose eigenvalues give the amplification factor
$$ g_j = 1 - 4r\sin^2\!\left(\frac{j\pi\Delta x}{2}\right) . $$
The most dangerous mode is the shortest wavelength the grid can carry — the sawtooth, alternating in sign from one point to the next — for which the sine equals one and $g = 1 - 4r$. Requiring $|g| \le 1$ for every mode again gives $r \le 1/2$.
The code below is the previous one with $r = 0.52$, a two per cent violation.
r = 0.52 # just above the limit
dt = r*dx**2/alpha
Nt = int(round(T/dt))
print('r = %.2f dt = %.2e number of time steps = %d' % (r, dt, Nt))
u = np.sin(np.pi*x/L)
u[0] = 0
u[-1] = 0
unew = np.zeros(Nx + 1)
print('\n step max error')
for n in range(Nt):
for i in range(1, Nx):
unew[i] = u[i] + r*(u[i+1] - 2*u[i] + u[i-1])
unew[0] = 0
unew[-1] = 0
u = unew.copy()
if (n + 1) % 100 == 0 or n + 1 == Nt:
print('%7d %10.2e' % (n + 1, np.max(np.abs(u - uexact(x, (n + 1)*dt)))))
# Plot
plt.plot(x, uexact(x, Nt*dt), 'k', linewidth = 2)
plt.plot(x, u, 'r', linewidth = 1.2)
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('u(x, t)', fontsize = 12)
plt.legend(['exact', 'r = 0.52'], fontsize = 11)
plt.show()
Two things are worth noticing. For the first four hundred steps the solution looks perfectly healthy — the error is no worse than in the stable run. Then it is destroyed, and what destroys it is the sawtooth, exactly the mode the analysis singled out. The instability was present from the very first step; it merely started from roundoff at $10^{-16}$ and needed a few hundred steps at a growth factor of $|1 - 4r| = 1.08$ to become visible. Raising $r$ to $0.55$ shortens that wait and the solution reaches $10^{18}$ instead.
The practical cost of the condition is severe, because it involves $\Delta x^2$. Halving the grid spacing to improve the spatial accuracy by a factor of four forces the time step down by a factor of four as well, so the total work grows by a factor of eight. Note also where the restriction comes from: the sawtooth, a mode that physically decays almost instantly and contributes nothing to the answer, dictates the step size for the whole computation. That is the definition of stiffness met in Section 4.4.2.
Chapter 4 had a cure for exactly this disease: evaluate the right-hand side at the new time level, which is backward Euler. Applied to the heat equation it gives the BTCS (backward time, central space) scheme,
$$ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha\left(\frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2}\right) , $$
or, after multiplying by $\Delta t$ and collecting the unknown new values on the left,
$$ -r\,u_{i-1}^{n+1} + (1 + 2r)\,u_i^{n+1} - r\,u_{i+1}^{n+1} = u_i^n , \qquad r = \frac{\alpha\,\Delta t}{\Delta x^2} . $$
Compare this with the explicit update of Section 6.2.1. There, the left-hand side contained one unknown and the right-hand side only known values. Here, three unknowns of the new level appear in every equation. Writing the equation out for $i = 1, 2, 3, \dots$ at a fixed time level $n$ shows the consequence:
$$ \begin{aligned} i = 1: \quad & -r\,u_0^{n+1} + (1+2r)\,u_1^{n+1} - r\,u_2^{n+1} = u_1^n, \\ i = 2: \quad & -r\,u_1^{n+1} + (1+2r)\,u_2^{n+1} - r\,u_3^{n+1} = u_2^n, \\ i = 3: \quad & -r\,u_2^{n+1} + (1+2r)\,u_3^{n+1} - r\,u_4^{n+1} = u_3^n, \\ & \;\;\vdots \\ i = N_x - 1: \quad & -r\,u_{N_x-2}^{n+1} + (1+2r)\,u_{N_x-1}^{n+1} - r\,u_{N_x}^{n+1} = u_{N_x-1}^n . \end{aligned} $$
The unknown $u_2^{n+1}$ appears in the first line, the second and the third; $u_3^{n+1}$ in the second, third and fourth; and so on. No line can be evaluated on its own, because each contains two unknowns that are only determined by its neighbours. The new values are coupled to each other, and the list of lines is not a loop but a linear system in the $N_x - 1$ unknowns $u_1^{n+1}, \dots, u_{N_x-1}^{n+1}$. This is the difference between explicit and implicit in one sentence: in the explicit scheme every line has one unknown and the inner loop just evaluates it; in the implicit scheme every line has three, and the inner loop is replaced by a solve.
Two of the lines contain a value that is not an unknown. The first contains $u_0^{n+1}$ and the last contains $u_{N_x}^{n+1}$, which are the boundary values at the new time level, given by the boundary conditions. As in Chapter 5, a known quantity is moved to the right-hand side, and the first and last lines become
$$ \begin{aligned} i = 1: \quad & (1+2r)\,u_1^{n+1} - r\,u_2^{n+1} = u_1^n + r\,u_0^{n+1}, \\ i = N_x - 1: \quad & -r\,u_{N_x-2}^{n+1} + (1+2r)\,u_{N_x-1}^{n+1} = u_{N_x-1}^n + r\,u_{N_x}^{n+1} . \end{aligned} $$
Now each line becomes one row of a matrix, with the coefficient of $u_j^{n+1}$ in column $j$: the entry $1+2r$ on the diagonal, and $-r$ on either side of it, except in the first and last rows where one of the $-r$'s has moved to the right with its boundary value. Each step therefore requires solving a linear system $A u^{n+1} = b$ ,
$$ \underbrace{ \begin{bmatrix} 1+2r & -r & & \\ -r & 1+2r & -r & \\ & \ddots & \ddots & \ddots \\ & & -r & 1+2r \end{bmatrix} }_{A} \underbrace{ \begin{bmatrix} u_1^{n+1} \\ u_2^{n+1} \\ \vdots \\ u_{N_x-1}^{n+1} \end{bmatrix} }_{u^{n+1}} = \underbrace{ \begin{bmatrix} u_1^n + r\,u_0^{n+1} \\ u_2^n \\ \vdots \\ u_{N_x-1}^n + r\,u_{N_x}^{n+1} \end{bmatrix} }_{b} . $$
This is precisely the kind of system built in Chapter 5 — tridiagonal, one row per interior point, with the boundary values entering the right-hand side. For the code, it is worth separating what depends on the time step from what does not.
The matrix $A$ contains only $r$. It does not change from step to step, so it is assembled once, outside the time loop: a diagonal of $1+2r$ and two off-diagonals of $-r$, of size $(N_x - 1) \times (N_x - 1)$.
The right-hand side $b$ is rebuilt at every step. Its interior entries are simply the current values $u_1^n, \dots, u_{N_x-1}^n$, that is, the interior of the current array; its first and last entries receive in addition $r$ times the boundary value at the new time level.
The step itself is one call to a linear solver, which returns all interior values of the new level at once. They are written into positions $1, \dots, N_x - 1$ of the new array, the boundary values into positions $0$ and $N_x$, and the new array becomes the current one, $n \leftarrow n + 1$. The inner loop of the explicit method has become "build $b$, solve $A u^{n+1} = b$"; the outer loop over time is unchanged.
The amplification factor of this scheme is $g_j = 1/\bigl(1 + 4r\sin^2(j\pi\Delta x/2)\bigr)$, which is smaller than one in magnitude for every $r$ and every mode. The implicit method is unconditionally stable: the time step may be chosen for accuracy alone.
r = 2.5 # five times the explicit limit
dt = r*dx**2/alpha
Nt = int(round(T/dt))
print('r = %.1f dt = %.2e number of time steps = %d' % (r, dt, Nt))
# the same tridiagonal assembly as in Chapter 5
A = np.zeros((Nx - 1, Nx - 1))
for i in range(Nx - 1):
A[i, i] = 1 + 2*r
if i > 0:
A[i, i - 1] = -r
if i < Nx - 2:
A[i, i + 1] = -r
u = np.sin(np.pi*x/L)
u[0] = 0
u[-1] = 0
plt.plot(x, u, 'k', linewidth = 2)
t = 0
for n in range(Nt):
b = u[1:-1].copy()
b[0] = b[0] + r*u[0] # the two boundary values
b[-1] = b[-1] + r*u[-1]
u[1:-1] = np.linalg.solve(A, b) # one linear solve per time step
t = t + dt
for tp in [0.01, 0.03, 0.05, 0.10]:
if abs(t - tp) < dt/2:
plt.plot(x, u, linewidth = 1.5)
plt.xlim([0, 1])
plt.xlabel('x', fontsize = 12)
plt.ylabel('u(x, t)', fontsize = 12)
plt.legend(['t = 0', 't = 0.01', 't = 0.03', 't = 0.05', 't = 0.10'], fontsize = 9, loc = 'upper right')
plt.show()
print('max error at t = %.2f %8.2e' % (t, np.max(np.abs(u - uexact(x, t)))))
# the whole history in one array: one row per time level, one column per grid point
U = np.zeros((Nt + 1, Nx + 1))
U[0, :] = np.sin(np.pi*x/L) # initial temperature (both ends already zero)
for n in range(Nt):
b = U[n, 1:-1].copy() # right-hand side from the current level
b[0] = b[0] + r*U[n+1, 0] # the two boundary values at the new level
b[-1] = b[-1] + r*U[n+1, -1]
U[n+1, 1:-1] = np.linalg.solve(A, b) # one linear solve per time step
tgrid = np.arange(Nt + 1)*dt
X, Tg = np.meshgrid(x, tgrid) # X[n, i] = x_i, Tg[n, i] = t_n, U[n, i] = u(x_i, t_n)
plt.figure(figsize = (9, 6))
ax = plt.axes(projection = '3d')
ax.plot_surface(X, Tg, U, cmap = 'viridis', rstride = 2, cstride = 1, linewidth = 0, alpha = 0.5)
for tp in [0, 0.01, 0.03, 0.05, 0.10]:
n = int(round(tp/dt)) # time level closest to tp
ax.plot(x, tp*np.ones(Nx + 1), U[n, :], linewidth = 2, color = 'r')
ax.set_xlabel('x', fontsize = 12)
ax.set_ylabel('t', fontsize = 12)
ax.set_zlabel('u(x, t)', fontsize = 12)
ax.view_init(elev = 30, azim = -45)
plt.show()
With $r = 2.5$, five times the explicit limit, the solution is perfectly well behaved: one hundred steps instead of five hundred and fifty-six, and no sign of a sawtooth. The price is accuracy rather than stability — the error is $1.9 \times 10^{-3}$ against $2.1 \times 10^{-4}$, because backward Euler is only first order in time and the step is five and a half times longer.
| explicit (FTCS) | implicit (BTCS) | |
|---|---|---|
| time derivative | forward difference | backward difference |
| space derivative evaluated at | $t_n$ | $t_{n+1}$ |
| work per step | one algebraic update | one linear solve $Au = b$ |
| stability | conditional, $r \le 1/2$ | unconditional, any $r > 0$ |
| cost | cheap per step, many steps | dearer per step, few steps |
That trade is what decides which method is used. If the accuracy wanted in time is comparable to the stability limit anyway, the explicit method is simpler and cheaper. If the stability limit is far below what accuracy requires — the usual situation on a fine grid — the implicit method wins, and the cost of a tridiagonal solve per step is repaid many times over.
Remark (other boundary conditions). An insulated end, $\partial u/\partial x = 0$ at $x = L$, is the Neumann condition of Section 5.2.3 and is handled in the same way: $u_{N_x}$ becomes an unknown, the system gains one row, and that row is a discretization of the boundary condition rather than of the PDE. The crude version $u_{N_x}^{n+1} = u_{N_x-1}^{n+1}$ is first-order accurate and costs an order at that end; the three-point one-sided formula preserves second order.
Remark (two space dimensions). On a rectangle, $u_{xx} + u_{yy}$ is discretized by applying the same three-point stencil in each direction, giving the five-point stencil $(u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 4u_{ij})/\Delta x^2$. Everything above carries over: the explicit method is still one line, its stability limit tightens to $r \le 1/4$, and the implicit method still requires a linear solve. What changes is the matrix. With $N$ points in each direction there are $N^2$ unknowns, and the matrix is no longer tridiagonal but sparse with five nonzero diagonals. Storing it as a full array, as was done here for clarity, becomes impossible almost immediately, and solving it efficiently — by sparse factorization or by iterative methods — is the subject on which large-scale scientific computing rests.
This chapter completes the pattern of the whole notebook. A derivative was replaced by a difference of neighbouring values in Chapter 2; that difference became a time-stepping rule in Chapter 4 and a matrix in Chapter 5; here it is both at once, and the balance between them — how far one may step in time given how finely one has divided space — is the central question of numerical simulation.
Key takeaways
Roots of nonlinear equations
- Roots of nonlinear equations are found iteratively. What a program can measure is the residual $f(\tilde r)$, a backward error: a small residual means that an exactly solvable nearby problem has been solved.
- Bisection needs only a sign change, $f(a) f(b) < 0$, and continuity. It cannot fail, and its error obeys the a-priori bound $|x_k - r| \le (b-a)/2^k$ — linear convergence with rate $1/2$, about one decimal digit per three evaluations.
- Fixed-point iteration $x_{k+1} = g(x_k)$ converges linearly with rate $|g'(p)|$ when this is less than $1$, and diverges otherwise; the rate is decided by how the equation was rearranged.
- Newton's method solves the tangent-line model at every step, $x_{k+1} = x_k - f(x_k)/f'(x_k)$. Because the rearrangement makes $g'(r) = 0$, the convergence to a simple root is quadratic and the number of correct digits roughly doubles per iteration — but only from a start close enough to the root.
- The secant method replaces $f'$ by a difference quotient through the last two iterates, costs one function evaluation per step, and converges with order $(1 + \sqrt 5)/2 \approx 1.618$, which usually makes it the cheapest of the three per digit gained.
- Guarantee and speed pull in opposite directions: bisection has the first and not the second, Newton the second and not the first. Practical software keeps a bracket and takes fast steps inside it.
Finite difference derivatives
- A finite difference is the definition of the derivative with the limit stopped short. One pair of Taylor expansions, one step to the right and one to the left, produces every formula of Chapter 2 and, at the same time, its error.
- Forward and backward differences are first-order accurate, with leading errors $\pm\tfrac{\Delta x}{2}f''$ — equal in size, opposite in sign. Plotted over an interval they appear shifted half a step to the left and to the right of the exact derivative.
- The central difference averages the two, cancels that term, and is second-order accurate at the same cost of two function evaluations. Symmetry is what buys the accuracy: a symmetric stencil cannot contain an odd power of $\Delta x$ in its error.
- Adding rather than subtracting the two expansions gives $f'' \approx (f_{k+1} - 2f_k + f_{k-1})/\Delta x^2$, the discrete Laplacian, also second-order accurate. Any order can be reached by using more points and choosing the coefficients so that the low-order error terms cancel.
- An order of accuracy is a slope on a log-log plot of error against $\Delta x$, and checking that slope is the standard test of a discretization.
- On a grid, use the central difference on the interior points and one-sided formulas at the two ends. A first-order end formula costs an order of accuracy there; a second-order one, built from three points, restores it.
Numerical integration
- Numerical integration replaces $\int_a^b f\,dx$ by a weighted sum of samples, $\sum_k w_k f(x_k)$. The methods differ only in the weights.
- The definition of the integral is itself the first method: with $\Delta x = (b-a)/N$ and $x_k = a + k\Delta x$, dropping the limit from the Riemann sum gives the left-sided rule $\Delta x\sum_{k=0}^{N-1} f_k$ and the right-sided rule $\Delta x\sum_{k=1}^{N} f_k$.
- Both rectangle rules are first-order accurate, with errors $\mp\tfrac{\Delta x}{2}(f(b) - f(a))$ — equal in size and opposite in sign. They differ from each other by exactly $(f(b) - f(a))\Delta x$.
- Averaging them cancels that term and gives the trapezoidal rule, $\tfrac{\Delta x}{2}\sum(f_k + f_{k+1})$, equivalently $\Delta x(\tfrac{f_0}{2} + f_1 + \cdots + f_{N-1} + \tfrac{f_N}{2})$. It is second-order accurate at no extra cost, with error $-\tfrac{\Delta x^2}{12}(f'(b) - f'(a))$ — and spectrally accurate for a smooth periodic integrand over a whole number of periods.
- Monte Carlo integration replaces the grid by random samples: with $x_j$ drawn uniformly from $\Omega$, the integral is the expected value $|\Omega|\,\mathbb{E}[f(X)]$ and is estimated by the sample mean $\frac{|\Omega|}{N}\sum_j f(x_j)$. The hit-or-miss variant does the same for an area, counting the fraction of points that land inside the region.
- Its error decays only as $O(N^{-1/2})$, far more slowly than a grid rule in one dimension, but that exponent is independent of dimension. A grid rule in $d$ dimensions has spacing $N^{-1/d}$, so the trapezoidal error becomes $O(N^{-2/d})$; the two rates cross at $d = 4$, and beyond it the grid rule loses badly. This is the curse of dimensionality, and escaping it is why sampling is used.
Numerical solutions to ODEs
- Any differential equation is put in the form $\dot{x} = f(x)$, $x(0) = x_0$, by adding derivatives as extra states. A numerical solver returns not a function but a trajectory: a list of states $x_0 \to x_1 \to \cdots \to x_N$ at the times $t_k = k\Delta t$.
- Forward Euler, $x_{k+1} = x_k + \Delta t f(x_k)$, evaluates the slope at the known point and is explicit. Backward Euler, $x_{k+1} = x_k + \Delta t f(x_{k+1})$, evaluates it at the unknown point and requires solving an equation each step. For a linear system they become $(I + \Delta t A)x_k$ and $(I - \Delta t A)^{-1}x_k$.
- The local error of one Euler step is $O(\Delta t^2)$, but $N \propto 1/\Delta t$ steps are needed, so the global error is $O(\Delta t)$: both methods are first order, and the measured error ratio under step halving is $2$. Forward and backward Euler are the rectangle rules of Chapter 3 in disguise, and averaging them gives the second-order Crank-Nicolson method.
- Stability is a separate question from accuracy. On $\dot y = \lambda y$ each step multiplies by a fixed factor, so with $z = \lambda\Delta t$ forward Euler is stable only for $|1 + z| \le 1$ — a disc of radius one centred at $-1$ — while backward Euler is stable for $|1 - z| \ge 1$, which includes the whole left half-plane.
- For a system, the same test reads $|\text{eig}(I + \Delta t A)| \le 1$. With real negative eigenvalues this caps the step at $\Delta t \le 2/|\lambda|_{\max}$, set by the fastest mode even when that mode has long since decayed. Systems with widely separated eigenvalues are stiff, and implicit methods exist to handle them.
- The same loop works unchanged for nonlinear $f$. On an undamped oscillator, whose eigenvalues are purely imaginary, forward Euler is unstable for every step size and steadily injects energy — a failure of the method, not of the step size.
Boundary value problems
- A boundary value problem gives conditions at the two ends of the interval, so there is no complete state to march forward from and the unknowns are all coupled. Replacing $u''$ by the three-point central stencil at every interior point turns $u'' = f$ into a linear system $Au = b$ with one equation per unknown.
- The matrix is the discrete Laplacian $(1, -2, 1)/\Delta x^2$ and is tridiagonal, because the stencil couples each point only to its neighbours. The known boundary values never enter the matrix; they move to the first and last entries of the right-hand side.
- The scheme inherits the order of its stencil: the error falls by a factor of $4$ each time the grid is halved. A Neumann condition adds one row, and a first-order boundary formula there would cost that second order.
- A time-stepping method applies the stencil sequentially; a boundary value method assembles it into a matrix and solves once.
Partial differential equations
- A PDE has several independent variables. Second-order equations are hyperbolic (waves), parabolic (diffusion) or elliptic (equilibrium) according to the sign of $B^2 - 4AC$; Chapter 5 was the elliptic case in one dimension, and the heat equation $u_t = \alpha u_{xx}$ is the parabolic model problem.
- Discretizing both variables gives $u_i^n$ on a grid in space and time. A forward difference in time and the central stencil in space give the explicit scheme $u_i^{n+1} = u_i^n + r(u_{i+1}^n - 2u_i^n + u_{i-1}^n)$ with mesh ratio $r = \alpha\Delta t/\Delta x^2$ — one line, no system to solve.
- The explicit scheme is stable only for $r \le 1/2$. Beyond that the centre weight $1 - 2r$ turns negative, the shortest grid mode is amplified by $|1 - 4r| > 1$, and the solution fills with a sawtooth — which grows from roundoff and so may stay invisible for hundreds of steps. Because the limit involves $\Delta x^2$, refining the grid by two multiplies the work by eight.
- Evaluating the space derivative at the new time level gives the implicit scheme, which requires solving a tridiagonal system — the one assembled in Chapter 5 — at every step, and is unconditionally stable. The choice between the two is cheap steps under a stability limit against dearer steps without one.
Looking ahead. One idea runs through all six chapters: replace the problem you cannot solve by a local model you can — a tangent or secant line for a root, a truncated Taylor series for a derivative, a constant or a straight line for an integral, one step of an integral for a differential equation. Each of the topics continues from where it stops here. For nonlinear equations the derivative becomes the Jacobian matrix and each Newton step is a linear system, and the same trade-off between a fast step and a safe one reappears in optimization. The finite difference method of Chapters 5 and 6 continues into two and three space dimensions, where the same stencil gives a large sparse matrix rather than a tridiagonal one, and solving it efficiently — rather than by the general-purpose elimination used here — is what makes large-scale simulation possible. Quadrature continues by interpolating at more points and integrating the interpolant, which is how the higher-order rules in reference books are generated. In the other direction it continues by sampling: the $O(N^{-1/2})$ rate that is indifferent to dimension is what makes expected values over high-dimensional distributions computable at all, and it underlies the stochastic estimates of gradients and losses on which machine learning rests. And for initial value problems, evaluating $f$ at several points inside a step buys much higher order for modest extra cost — the Runge-Kutta family, whose fourth-order member is the engine inside solve_ivp's default RK45, with adaptive step size control on top.