This is a Jupyter notebook using Python. You can install Jupyter locally to edit and interact with this notebook.
Consider the boundary value problem
$$ \begin{gather} -\frac{d^2 u}{dx^2} = f(x) \quad x \in \Omega = (-1,1) \\ u(-1) = a \quad \frac{du}{dx}(1) = b . \end{gather} $$$f(x)$ is the "forcing" term and we have a Dirichlet boundary condition at the left end of the domain and a Neumann condition on the right end. We need to choose
How can we compute $\frac{du}{dx}(x_i)$ using $u_i$ and $u_{i+1}$? How accurate is this approximation?
Without loss of generality, we may assume $x_i = 0$ by shifting our function. For notational convenience, we will also define $h = x_{i+1} - x_i$. To determine the order of accuracy, we represent the function $u(x)$ by its Taylor series $$ u(x) = u(0) + u'(0)x + u''(0)x^2/2! + O(x^3)$$ and substitute into the differencing formula $$ \begin{split} u'(0) \approx \frac{u(h) - u(0)}{h} = h^{-1} \Big( u(0) + u'(0) h + u''(0)h^2/2 + O(h^3) - u(0) \Big) \\ = u'(0) + u''(0)h/2 + O(h^2) . \end{split}$$ Evidently the error in this approximation is $u''(0)h/2 + O(h^2)$. We say this method is first order accurate.
In [ ]:
In [1]:
%matplotlib notebook
import numpy
from matplotlib import pyplot
pyplot.style.use('ggplot')
n = 20
h = 2/(n-1)
x = numpy.linspace(-1,1,n)
u = numpy.sin(x)
pyplot.figure()
pyplot.plot(x[:-1], (u[1:] - u[:-1]) / h, 'o')
pyplot.plot(x[:-1], numpy.cos(x[:-1]))
Out[1]:
This "has the right shape", but the numerical approximation is "shifted" to the left of the analytic derivative. If we differenced to the left, it would be shifted the other way.
Here we try
$$ \begin{split} u'(0) \approx \frac{u(h) - u(-h)}{2h} \\ = (2h)^{-1} \Big( \big[ u(0) + u'(0)h + u''(0)h^2/2 + u'''(0)h^3/6 + O(h^4) \Big] - \big[ u(0) - u'(0)h + u''(0)h^2/2 - u'''(0)h^3/6 + O(h^4) \big] \Big) \\ = u'(0) + u'''(0)h^2/6 + O(h^3) \end{split} $$
In [2]:
def diff1l(x, u):
return x[:-1], (u[1:] - u[:-1]) / (x[1:] - x[:-1])
def diff1r(x, u):
return x[1:], (u[1:] - u[:-1]) / (x[1:] - x[:-1])
def diff1c(x, u):
return x[1:-1], (u[2:] - u[:-2]) / (x[2:] - x[:-2])
pyplot.figure()
for diff in (diff1l, diff1r, diff1c):
xx, yy = diff(x, u)
pyplot.plot(xx, yy, 'o', label=diff.__name__)
pyplot.plot(x, numpy.cos(x), label='exact')
pyplot.legend(loc='upper right')
Out[2]:
In [3]:
grids = 2**numpy.arange(3,10)
def grid_refinement_error(f, fp, diff):
error = []
for n in grids:
x = numpy.linspace(-1, 1, n)
xx, yy = diff(x, f(x))
error.append(numpy.linalg.norm(yy - fp(xx), numpy.inf))
return grids, error
pyplot.figure()
for diff in (diff1l, diff1r, diff1c):
ns, error = grid_refinement_error(numpy.sin, numpy.cos, diff)
pyplot.loglog(ns, error, 'o', label=diff.__name__)
pyplot.loglog(grids, grids**(-1.), label='$1/n$')
pyplot.loglog(grids, grids**(-2.), label='$1/n^2$')
pyplot.xlabel('Resolution $n$')
pyplot.ylabel('Derivative error')
pyplot.legend(loc='upper right')
Out[3]:
In [4]:
x = numpy.linspace(-1, 1, 9)
xfine = numpy.linspace(-1, 1, 200)
def f_rough(x):
return numpy.cos(4*numpy.pi*x)
def fp_rough(x):
return -4*numpy.pi * numpy.sin(4*numpy.pi*x)
pyplot.figure()
pyplot.plot(x, f_rough(x), 'o-')
pyplot.plot(xfine, f_rough(xfine), '-')
Out[4]:
In [5]:
pyplot.figure()
for diff in (diff1r, diff1l, diff1c):
xx, yy = diff(x, f_rough(x))
pyplot.plot(xx, yy, 'o', label=diff.__name__)
pyplot.plot(xfine, fp_rough(xfine), label='exact')
pyplot.legend(loc='upper right')
Out[5]:
We have this function $f_{\text{rough}}(x)$ that is rough at the grid scale and has derivative identically zero according to the centered difference formula. Except perhaps for boundary conditions (which we'll consider later), given any solution $u(x)$ to a differential equation discretized using the centered difference approximation, $$\tilde u(x) = u(x) + f_{\text{rough}}(x)$$ would also be a solution. This non-uniqueness is a disaster for numerical algorithms and in most cases, the centered difference formula is not used directly to compute first derivatives.
When assessing a discretization we require both accuracy and stability.
We will need at least three grid points to compute a second derivative.
Again, there is more than one possible approach.
Why should we choose one over the other? Can we understand this in terms of accuracy and stability?
In [6]:
def diff2a(x, u):
xx, yy = diff1c(x, u)
return diff1c(xx, yy)
def diff2b(x, u):
xx, yy = diff1l(x, u)
return diff1r(xx, yy)
x = numpy.linspace(-1, 1, 11)
u = numpy.cos(x)
pyplot.figure()
for diff2 in (diff2a, diff2b):
xx, yy = diff2(x, u)
pyplot.plot(xx, yy, 'o', label=diff2.__name__)
pyplot.plot(x, -numpy.cos(x), label='exact')
pyplot.legend(loc='upper right')
Out[6]:
In [7]:
def grid_refinement_error2(f, fpp, diff):
error = []
for n in grids:
x = numpy.linspace(-1, 1, n)
xx, yy = diff(x, f(x))
error.append(numpy.linalg.norm(yy - fpp(xx), numpy.inf))
return grids, error
pyplot.figure()
for diff2 in (diff2a, diff2b):
ns, error = grid_refinement_error2(numpy.cos, lambda x: -numpy.cos(x), diff2)
pyplot.loglog(ns, error, 'o', label=diff.__name__)
pyplot.loglog(grids, (grids-1)**(-1.), label='$h$')
pyplot.loglog(grids, (grids-1)**(-2.), label='$h^2$')
pyplot.loglog(grids, .25*(grids-1)**(-2.), label='$h^2/4$')
pyplot.xlabel('Resolution $n$')
pyplot.ylabel('Second derivative error')
pyplot.legend(loc='upper right')
Out[7]:
diff2a method is more accurate by a factor of 4diff2b method cannot compute the derivative at points next to the boundarySo far we have written functions of the form diff(x, u) that compute derivatives. These functions happen to have been linear in u. We should be able to write differentiation as a matrix $D$ such that
where $x$ is the vector of $n$ discrete points, thus $u(x)$ is also a vector of length $n$.
hw0 inside the repository. Add your source file(s) to that directory.diffmat(x) that returns a matrix $D$ that computes first derivatives.diff2mat(x) that returns a matrix $D_2$ that computes second derivatives.README.md in the hw0 directory and summarize your results (one paragraph or a few bullet items is fine).x are monotonically increasing.Our boundary value problem states a differential equation to be satisfied in the interior of the domain, combined with Dirichlet conditions at the left endpoint and Neumann at the right endpoint.
The left endpoint in our example BVP has a Dirichlet boundary condition, $$u(-1) = a . $$ With finite difference methods, we have an explicit degree of freedom $u_0 = u(x_0 = -1)$ at that endpoint. When building a matrix system for the BVP, we can implement this boundary condition by modifying the first row of the matrix, $$ \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ \\ & & A_{1:,:} & & \\ \\ \end{bmatrix} \begin{bmatrix} u_0 \\ \\ u_{1:} \\ \\ \end{bmatrix} = \begin{bmatrix} a \\ \\ f_{1:} \\ \\ \end{bmatrix} . $$
We can eliminate $u_0$ and create a reduced system for $u_{1:}$. We will describe this for a more general $2\times 2$ block system of the form $$ \begin{bmatrix} I & 0 \\ A_{10} & A_{11} \end{bmatrix} \begin{bmatrix} u_0 \\ u_1 \end{bmatrix} = \begin{bmatrix} f_0 \\ f_1 \end{bmatrix} .$$ We can rearrange as $$ A_{11} u_1 = f_1 - A_{10} f_0 $$ which is symmetric if $A_{11}$ is. This is sometimes called "lifting" and is often done implicitly in the mathematics literature. It is convenient for linear solvers and eigenvalue solvers, but inconvenient for IO and postprocessing, as well as some nonlinear problems. For this reason, it may be preferable to write $$ \begin{bmatrix} I & 0 \\ 0 & A_{11} \end{bmatrix} \begin{bmatrix} u_0 \\ u_1 \end{bmatrix} = \begin{bmatrix} f_0 \\ f_1 - A_{10} f_0 \end{bmatrix} $$ which is symmetric and entirely decouples the degrees of freedom associated with the boundary. This method turns out to be relatively elegant for nonlinear solvers.
It is sometimes useful to scale the identity by some scalar related to the norm of $A_{11}$.
The right endpoint in our example BVP has a Neumann boundary condition, $$ \frac{du}{dx}(1) = b . $$ With finite difference methods, there are generally two ways to derive an equation for the value $u_n = u(x_n = 1)$. The first is to use a one-sided difference formula as in $$ \frac{u_n - u_{n-1}}{h} = b . $$ This requires an extra discretization choice and it may not be of the same order of accuracy as the interior discretization and may destroy symmetry.
An alternative is to temporarily introduce a ghost value $u_{n+1} = u(x_{n+1} = 1 + h)$ (possibly more) and define it to be a reflection of the values from inside the domain. In the case $b=0$, this reflection is $u_{n+i} = u_{n-i}$. More generally, we can bias the reflection by the slope $b$, $$ u_{n+i} = u_{n-i} + 2b(x_n - x_{n-i}) . $$ After this definition of ghost values, we apply the interior discretization at the boundary. For our reference equation, we would write
$$ \frac{-u_{n-1} + 2 u_n - u_{n+1}}{h^2} = f(x_n) $$which simplifies to $$ \frac{u_n - u_{n-1}}{h^2} = f(x_n)/2 + b/h $$ after dividing by 2 and moving the boundary term to the right hand side.
In [44]:
def laplacian(n, rhsfunc, a, b):
x = numpy.linspace(-1, 1, n+1)
h = 2 / n
rhs = rhsfunc(x)
e = numpy.ones(n)
L = (2 * numpy.eye(n+1) - numpy.diag(e, 1) - numpy.diag(e, -1)) / h**2
# Dirichlet condition
L[0,:] = numpy.eye(1,n+1)
rhs[0] = a
rhs[1:] -= L[1:,0] * a
L[1:,0] = 0
# Neumann condition
L[-1,-1] /= 2
rhs[-1] = b / h + rhs[-1] / 2
return x, L, rhs
laplacian(5, lambda x: 0*x+1, 1, -1)
Out[44]:
In [24]:
x, L, rhs = laplacian(20, lambda x: 0*x+1, 1, -.5)
u = numpy.linalg.solve(L, rhs)
pyplot.figure()
pyplot.plot(x, u)
Out[24]:
A continuous Green's function is the solution of a differential equation where the right hand side is a Dirac delta function. We can compute discrete Green's functions by solving against columns of the identity.
In [10]:
Linv = numpy.linalg.inv(L)
pyplot.figure()
pyplot.plot(x, Linv[:,[5,10,15]])
Out[10]:
In [35]:
lam, U = numpy.linalg.eigh(L)
pyplot.figure()
pyplot.semilogy(lam, 'o')
Out[35]:
In [34]:
x, L, rhs = laplacian(80, lambda x: 0*x+1, 1, -.5)
lam, U = numpy.linalg.eigh(L)
pyplot.figure()
for i in (0, 1, 2, 3):
pyplot.plot(x, U[:,i], label='$\lambda_{%d} = %5f$'%(i, lam[i]))
pyplot.legend(loc='upper right')
Out[34]:
Above, we see evidence of stability in that all the eigenvalues are positive and the eigenvectors associated with the smallest eigenvalues are low frequency (ignoring the mode corresponding to the Dirichlet boundary). Specifically, there are no "noisy" functions with eigenvalues close to zero.
In [13]:
pyplot.figure()
for i in (-3, -2, -1):
pyplot.plot(x, U[:,i], label='$\lambda_{%d} = %5f$'%(i, lam[i]))
pyplot.legend(loc='upper right')
Out[13]:
In [45]:
def uexact(x):
return numpy.tanh(3*x)
def duexact(x):
return 3*numpy.cosh(3*x)**(-2)
def negdduexact(x):
return 3**2 * 2 * numpy.tanh(3*x) * numpy.cosh(3*x)**(-2)
x, L, f = laplacian(20, negdduexact, uexact(-1), duexact(1))
u = numpy.linalg.solve(L, f)
pyplot.figure()
pyplot.plot(x, uexact(x), label='exact')
pyplot.plot(x, u, 'o', label='computed')
pyplot.legend(loc='upper left')
Out[45]:
In [46]:
def mms_error(n, discretize):
x, L, f = discretize(n, negdduexact, uexact(-1), duexact(1))
u = numpy.linalg.solve(L, f)
return numpy.linalg.norm(u - uexact(x), numpy.inf)
ns = 2**numpy.arange(3,8)
errors = [mms_error(n, laplacian) for n in ns]
pyplot.figure()
pyplot.loglog(2/ns, errors, 'o', label='numerical')
pyplot.loglog(2/ns, (2/ns), label='$h$')
pyplot.loglog(2/ns, (2/ns)**2, label='$h^2$')
pyplot.legend(loc='upper left')
Out[46]:
We have been using the discretization
$$ -u''(x_i) \approx \frac{-u_{i-1} + 2 u_i - u_{i+1}}{h^2} $$which is sometimes summarized using the stencil
$$ h^{-2} \begin{bmatrix} -1 & 2 & -1 \end{bmatrix} . $$There is also the discretization arising from applying centered first derivatives twice,
$$ -u''(x_i) \approx \frac{-u_{i-2} + 2 u_i - u{i+2}}{4 h^2} $$which corresponds to the stencil
$$ (2h)^{-2} \begin{bmatrix} -1 & 0 & 2 & 0 & -1 \end{bmatrix} . $$We'll reuse the boundary condition handling and keep the first stencil at the points adjacent to the boundaries, $x_1$ and $x_{n-1}$, but apply this discretization for true interior points.
In [47]:
def laplacian2(n, rhsfunc, a, b):
x, L, rhs = laplacian(n, rhsfunc, a, b)
h = 2 / n
L[2:-2,:] = (2 * numpy.eye(n-3, n+1, k=2)
- numpy.eye(n-3, n+1, k=0)
- numpy.eye(n-3, n+1, k=4)) / (2*h)**2
return x, L, rhs
laplacian2(6, numpy.cos, 1, 0)
Out[47]:
This matrix is not symmetric due to the discretization used near the boundaries.
In [48]:
x, L, f = laplacian2(20, negdduexact, uexact(-1), duexact(1))
u = numpy.linalg.solve(L, f)
pyplot.figure()
pyplot.plot(x, uexact(x), label='exact')
pyplot.plot(x, u, 'o', label='computed')
pyplot.legend(loc='upper left')
Out[48]:
In [49]:
errors = [mms_error(n, laplacian) for n in ns]
errors2 = [mms_error(n, laplacian2) for n in ns]
pyplot.figure()
pyplot.loglog(2/ns, errors, 'o', label='numerical')
pyplot.loglog(2/ns, errors2, 's', label='numerical2')
pyplot.loglog(2/ns, (2/ns), label='$h$')
pyplot.loglog(2/ns, (2/ns)**2, label='$h^2$')
pyplot.legend(loc='upper left')
Out[49]:
In [19]:
n = 40
x, L, _ = laplacian(n, lambda x: 0*x+1, 1, -1)
_, L2, _ = laplacian2(n, lambda x: 0*x+1, 1, -1)
lam, U = numpy.linalg.eigh(L)
lam2, U2 = numpy.linalg.eig(L2) # Need to use nonsymmetric eigensolver
idx = numpy.argsort(lam2)
lam2 = lam2[idx]
U2 = U2[:,idx]
pyplot.figure()
pyplot.plot(lam, 'o', label='original')
pyplot.plot(lam2, 'o', label='version 2')
pyplot.legend(loc='lower right')
Out[19]:
In [20]:
pyplot.figure()
for i in (0, 1, 2, 3):
pyplot.plot(x, U[:,i], label='$\lambda_{%d} = %5f$'%(i, lam[i]))
pyplot.legend(loc='upper right')
pyplot.title('Original')
pyplot.figure()
for i in (0, 1, 2, 3):
pyplot.plot(x, U2[:,i], label='$\lambda_{%d} = %5f$'%(i, lam2[i]))
pyplot.legend(loc='upper right')
pyplot.title('Version 2')
Out[20]:
laplacian2laplacian2 converges at the same order of accuracy as laplacianlaplacian2 is less accurate, especially on coarse/under-resolved gridslaplacian2 (with this boundary condition implementation) is nonsingularlaplacian2 is a weakly unstable discretization, even away from boundaries