26 P3: Hints

Problem Class 3: Finite Differences

Published

April 6, 2026

Abstract

Exercises based on content from Lectures 12-15

using Plots, LinearAlgebra, LaTeXStrings

Exercise 26.1 So far, we have considered homogeneous boundary conditions u(0) = u(1) = 0. How would you approximate the problem

\begin{align} &-u'' + p u' + q u = f \qquad \text{on }(0,1) \nonumber\\ &u(0) = \alpha, \qquad u(1) = \beta \end{align}

with finite differences?

Again, suppose we have a uniform grid x_j := j h for j=0,\dots,n and h := \frac{1}{n}.

For simplicity, we set p=0. Before, we approximated -u''(x_j) with the finite difference

\begin{align} (-Lu)_j = \frac{-u_{j-1} + 2 u_j - u_{j+1}}{h^2} \end{align}

and set u_0 = u_n = 0 to enforce the boundary condition. Now, we simply set u_0 := \alpha and u_n := \beta which leads to the matrix equation

\begin{align} (-L + q) u = \begin{pmatrix} f(x_1) + \frac{\alpha}{h^2}\\ f(x_2) \\ \vdots \\ f(x_{n-2})\\f(x_{n-1})\\f(x_n) + \frac{\beta}{h^2} \end{pmatrix}. \end{align}

Exercise 26.2 So far, we have considered ADR

\begin{align} -u'' + p u' + q u = f \end{align}

with constant coefficients p,q. How would you write down a finite difference approximation to this equation if p and q are functions of x?

Before, we approximated p u'(x_j) \approx p (Du)_j (for some differentiation matrix D). Now, one may simply approximate p(x_j) u'(x_j) with p(x_j) (Du)_j. Moreover, we would replace q I with the diagonal matrix formed by evaluating q on the interior grid points. As before, this defines a linear system of equations that may be solved for the interior grid points.

Exercise 26.3 (FNC example 10.4.2) Following on from the previous two questions, we try and solve the problem

\begin{align} &-u''(x) + a \cos(ax) u'(x) - a^2\sin(ax) u(x) = 0\nonumber\\ &u(0) = 1 \nonumber\\ &u(\tfrac{3\pi}{2a}) = e^{-1}. \end{align}

Show that the exact solution is u(x) = e^{\sin ax}. Read through and try to understand the following code:

"""
    adr_1d(n, p, q, f; a=0,b=1, α=0,β=0)

Approximate the solution to 
    -u′′ + p u′ + q u = 0 on (a,b)
    u(a) = α, u(b) = β
(linear Advection-Diffusion-Reaction equation with Dirichlet boundary conditions)
using central finite differences
"""
function adr_1d(n, p, q, f; a=0.,b=1., α=0.,β=0.)
    # mesh size
    h = (b-a)/n 
    # interior grid points 
    x = [ a + j *(b-a)/n for j ∈ 1:n-1 ] 

    # Finite Difference Coefficients
    # (using a second order central difference) 
    d = (2/h^2) * ones(n-1) + q.(x)
    d₋  = (-1/h^2)*ones(n-2) - p.(x[2:end])/(2h)
    d₊  = (-1/h^2)*ones(n-2) + p.(x[1:end-1])/(2h)
    Lh = Tridiagonal(d₋, d, d₊)
    F = f.(x)
    F[1] += (1/h^2 + p(x[1])/(2h))*α
    F[end] += (1/h^2 - p(x[end])/(2h))*β
    U = Lh \ F

    # add the boundary conditions
    U = [α, U..., β]
    x = [a, x..., b]

    return x, U
end

p(x) = cos(x)
q(x) = -sin(x)
f(x) = 0.    

x, u = adr_1d(50, p, q, f; a=0.,b=3π/2, α=1.,β=exp(-1))
P1 = scatter(x, u, title=L"$-u′′ + \cos(x) u′ - \sin(x) u = 0$", 
     label="finite difference approximation", markersize=2)
u_e(x) = exp(sin(x))
plot!( u_e, label="exact solution" )

Firstly, notice that u(x) = e^{\sin ax} satisfies u(0) = e^0 = 1 and u(\tfrac{3\pi}{2a}) = e^{sin \tfrac{3\pi}{2}} = 1/e. Moreover,

\begin{align} u'(x) &= a \cos(ax) e^{\sin ax} \nonumber\\ &= a \cos(ax) u(x)\qquad \text{and}\\ u''(x) &= -a^2\sin(ax) u(x) + a\cos(ax) u'(x) \nonumber\\ &= a^2\left[ - \sin(ax) + \cos^2(ax) \right] u(x) \end{align}

As a result,

\begin{align} &-u''(x) + a \cos(ax) u'(x)\nonumber\\ % &= a^2\left[ \sin(ax) - \cos^2(ax) \right] u(x) % + a^2 \cos^2(ax) u(x) \nonumber\\ % &= a^2 \sin(ax) u(x) \end{align}

as required.

Exercises from the end of lecture 15

Exercise 26.4 (i) Show that \hat{U}_0 = \frac{1}{\sqrt{n}} \sum_{j=0}^{n-1} U_j.
(ii) Write -L_{\rm per} U = F in terms of \hat{U}_k and \hat{F}_k,
(iii) Hence show that, if F has zero mean \sum_{j=0}^{n-1} F_j = 0, then there exists U solving -L_{\rm per} U = F.

Recall that \bm e_k defined by (\bm e_k)_j := \frac{1}{\sqrt{n}} e^{\frac{2\pi \mathrm{i}jk}{n}} for k=0,\dots,n-1 forms an orthonormal basis of \mathbb C^{n} and we defined \hat{U}_k as the coefficients of U in this basis: U_j = \sum_{k=0}^{n-1} \hat{U}_k (\bm e_k)_j. Therefore, using the orthonormality, we have

\begin{align} \hat{U}_k = \sum_{j=0}^{n-1} U_j \overline{(\bm e_k)_j}. \end{align}

(i) Chosing k=0 gives \hat{U}_0 = \sum_{j=0}^{n-1} U_j \overline{(\bm e_0)_j} = \frac{1}{\sqrt{n}}\sum_{j=0}^n U_j.

(ii) Notice that

\begin{align} (-L_{\rm per} U)_j &= \sum_{k=0}^{n-1} \hat{U}_k (-L_{\rm per}\bm e_k)_j \nonumber\\ % &= \sum_{k=0}^{n-1} \hat{U}_k \frac{-(\bm e_k)_{j-1} + 2(\bm e_k)_j - (\bm e_k)_{j+1}}{h^2} \nonumber\\ % &= \sum_{k=0}^{n-1} \hat{U}_k \frac {2 - \big( e^{-\frac{2\pi\mathrm{i}k}{n}} + e^{\frac{2\pi\mathrm{i}k}{n}} \big)} {h^2} (\bm e_k)_j \nonumber\\ % &= \sum_{k=0}^{n-1} \hat{U}_k \frac {2 - 2\cos\frac{2\pi k}{n}} {h^2} (\bm e_k)_j. \end{align}

The right hand side is F_j = \sum_{k=0}^{n-1} \hat{F}_k (\bm e_k)_j.

(iii) By the previous part, we want to solve

\begin{align} \hat{U}_k \frac {2 - 2\cos\frac{2\pi k}{n}} {h^2} % &= \hat{F}_k. \end{align}

for k=0,\dots,n-1. In particular, when k=0, we need \hat{F}_0 = \hat{U}_0 \frac {2 - 2\cos 0} {h^2} = 0. By the first part, we have \hat{F}_0 = \frac{1}{\sqrt{n}}\sum_{j=0}^{n-1} F_j. Moveover, 2 - 2\cos\frac{2\pi k}{n} \not= 0 for all k=1,\dots,n-1 and so we may solve \hat{U}_k \frac {2 - 2\cos\frac{2\pi k}{n}} {h^2} % = \hat{F}_k for all k=1,\dots,n-1 (by dividing by \lambda_k := \frac {2 - 2\cos\frac{2\pi k}{n}} {h^2} \not=0).

If \sum_{j=0}^{n-1} F_j = 0, the general solution to -L_{\rm per} U = F is given by U = C \bm e_0 + \sum_{k=1}^{n-1} \frac{\hat{F}_k}{\lambda_k} \bm e_k where C is an arbitrary constant (equal to the mean of U).

Exercise 26.5 Read the introduction to the book

“Brigham, E. O. (1988). The Fast Fourier Transform and Its Applications.”

Other Links