Residual#
We need to have a way to tell when to stop smoothing.
If we know the analytic solution, then we can just compare to that, and keep iterating until the error is small, but that kinda defeats the purpose. Instead, we can measure how well we satisfy the discrete equation—this is called the residual.
We still need something to compare to, so we define the source norm, \(\| f \|\), and we will pick a tolerance \(\epsilon\) and iterate until:
For the special case of a homogeneous source (\(f = 0\)), then we will iterate until
We will use the L2 norm:
Let’s update our grid class to include a norm and residual function:
import numpy as np
import matplotlib.pyplot as plt
This is the same grid class as before
class Grid:
def __init__(self, nx, ng=1, xmin=0, xmax=1,
bc_left_type="dirichlet", bc_left_val=0.0,
bc_right_type="dirichlet", bc_right_val=0.0):
self.xmin = xmin
self.xmax = xmax
self.ng = ng
self.nx = nx
self.bc_left_type = bc_left_type
self.bc_left_val = bc_left_val
self.bc_right_type = bc_right_type
self.bc_right_val = bc_right_val
# python is zero-based. Make easy integers to know where the
# real data lives
self.ilo = ng
self.ihi = ng+nx-1
# physical coords -- cell-centered
self.dx = (xmax - xmin)/(nx)
self.x = xmin + (np.arange(nx+2*ng)-ng+0.5)*self.dx
# storage for the solution
self.phi = self.scratch_array()
self.f = self.scratch_array()
def scratch_array(self):
"""return a scratch array dimensioned for our grid """
return np.zeros((self.nx+2*self.ng), dtype=np.float64)
def norm(self, e):
"""compute the L2 norm of e that lives on our grid"""
return np.sqrt(self.dx * np.sum(e[self.ilo:self.ihi+1]**2))
def fill_bcs(self):
"""fill the boundary conditions on phi"""
# we only deal with a single ghost cell here
# left
if self.bc_left_type.lower() == "dirichlet":
self.phi[self.ilo-1] = 2 * self.bc_left_val - self.phi[self.ilo]
elif self.bc_left_type.lower() == "neumann":
self.phi[self.ilo-1] = self.phi[self.ilo] - self.dx * self.bc_left_val
else:
raise ValueError("invalid bc_left_type")
# right
if self.bc_right_type.lower() == "dirichlet":
self.phi[self.ihi+1] = 2 * self.bc_right_val - self.phi[self.ihi]
elif self.bc_right_type.lower() == "neumann":
self.phi[self.ihi+1] = self.phi[self.ihi] - self.dx * self.bc_right_val
else:
raise ValueError("invalid bc_right_type")
Now we’ll write our residual routine
def residual_norm(g):
"""compute the residual norm"""
r = g.scratch_array()
r[g.ilo:g.ihi+1] = g.f[g.ilo:g.ihi+1] - (g.phi[g.ilo+1:g.ihi+2] -
2 * g.phi[g.ilo:g.ihi+1] +
g.phi[g.ilo-1:g.ihi]) / g.dx**2
return g.norm(r)
Now we’ll write a relaxation function that does smoothing until either a maximum number of iterations is taken or we reach a desired tolerance. If the tolerance is set to None, then the routine will take the full amount of iterations.
class TooManyIterations(Exception):
pass
def relax(g, tol=1.e-8, max_iters=200000, analytic=None):
iter = 0
fnorm = g.norm(g.f)
# if the RHS is f = 0, then set the norm to 1
# this ensures that we then check for residual < eps
if fnorm == 0.0:
fnorm = 1
g.fill_bcs()
r = residual_norm(g)
res_norm = []
true_norm = []
while iter < max_iters and (tol is None or r > tol * fnorm):
g.phi[g.ilo:g.ihi+1:2] = 0.5 * (-g.dx * g.dx * g.f[g.ilo:g.ihi+1:2] +
g.phi[g.ilo+1:g.ihi+2:2] + g.phi[g.ilo-1:g.ihi:2])
g.fill_bcs()
g.phi[g.ilo+1:g.ihi+1:2] = 0.5 * (-g.dx * g.dx * g.f[g.ilo+1:g.ihi+1:2] +
g.phi[g.ilo+2:g.ihi+2:2] + g.phi[g.ilo:g.ihi:2])
g.fill_bcs()
r = residual_norm(g)
res_norm.append(r / fnorm)
if analytic is not None:
true_norm.append(g.norm(g.phi - analytic(g.x)))
iter += 1
if tol is not None and iter >= max_iters:
raise TooManyIterations(f"too many iterations, niter = {iter}")
return res_norm, true_norm, iter
Residual vs. Truncation Error#
Now we’ll look at how the residual error compares to the truncation error of our discretization. We’ll take a fixed number of iterations for our same model problem,
on \([0, 1]\) with homogeneous Dirichlet BCs.
def analytic(x):
return -np.sin(x) + x * np.sin(1.0)
def f(x):
return np.sin(x)
from matplotlib.lines import Line2D
fig, ax = plt.subplots()
for i, nx in enumerate([16, 32, 64]):
g = Grid(nx)
g.f[:] = f(g.x)
res_norm, true_norm, _ = relax(g, tol=None, max_iters=20000, analytic=analytic)
n = np.arange(len(res_norm)) + 1
ax.loglog(n, true_norm, label=f"{nx}", color=f"C{i}")
ax.loglog(n, res_norm, ls=":", color=f"C{i}")
print(f"nx = {nx}, true error = {true_norm[-1]}")
# manually add line types to the legend
handles, labels = ax.get_legend_handles_labels()
line1 = Line2D([0], [0], ls="-", color="k", label="truncation error")
line2 = Line2D([0], [0], ls=":", color="k", label="residual")
handles.extend([line1, line2])
ax.legend(handles=handles, ncol=2)
ax.set_xlabel("number of smoothing iterations")
ax.grid(ls=":", color="0.5")
nx = 16, true error = 0.0002489743563076079
nx = 32, true error = 6.22480653738797e-05
nx = 64, true error = 1.5562295551083177e-05
Look at what this shows us:
The truncation error (solid line) stalls at a much higher value than the residual (dotted).
The truncation error converges second order with the number of zones (look at the numbers printed during the run)
As we increase the number of zones, we need more iterations until the residual drops to machine roundoff
The residual error eventually reaches roundoff—this indicates that we satisfy the discrete equation “exactly”
Note that the residual is telling us how well we satisfy the discrete form of Poisson’s equation, but it is not telling us anything about whether the discrete approximation we are using for Poisson’s equation is a good approximation. That is our job as a computational scientist to assess.
We are using a second-order accurate discretization. If needed, we could use a higher-order accurate method, although that would be more computationally expensive.
In-class exercise:
Solve the problem
on \([0, 1]\) with homogeneous Neumann boundary conditions.
C++ implementation#
Here’s a C++ implementation that supports Dirichlet boundary conditions.
This can be compiled as:
g++ -I. -o poisson poisson.cpp
and it will output a text file, poisson.txt with 3 columns of information: \(x\), \(\phi\), and \(f\).