# 2021-08-30 High order and boundaries¶

## Last time¶

• Stability

• Derivatives as matrices

• Boundary conditions

• Discrete eigenvalues/eigenvectors

## Today¶

• High order on arbitrary grids

• Method of manufactured solutions

• More satisfying approach to boundary conditions

using Plots


# A few methods on grids¶

diff1l(x, u) = x[2:end],   (u[2:end] - u[1:end-1]) ./ (x[2:end] - x[1:end-1])
diff1r(x, u) = x[1:end-1], (u[2:end] - u[1:end-1]) ./ (x[2:end] - x[1:end-1])
diff1c(x, u) = x[2:end-1], (u[3:end] - u[1:end-2]) ./ (x[3:end] - x[1:end-2])
difflist = [diff1l, diff1r, diff1c]

n = 20
h = 2 / (n - 1)
x = LinRange(-3, 3, n)
u = sin.(x)
fig = plot(cos, xlims=(-3, 3))
for d in difflist
xx, yy = d(x, u)
plot!(fig, xx, yy, marker=:circle, label=d)
end

fig # Measuring error on grids¶

using LinearAlgebra

grids = 2 .^ (2:10)
hs = 1 ./ grids
function refinement_error(f, fprime, d)
error = []
for n in grids
x = LinRange(-3, 3, n)
xx, yy = d(x, f.(x))
push!(error, norm(yy - fprime.(xx), Inf))
end
error
end

refinement_error (generic function with 1 method)

fig = plot(xscale=:log10, yscale=:log10)
for d in difflist
error = refinement_error(sin, cos, d)
plot!(fig, hs, error, marker=:circle, label=d)
end
plot!(fig, hs, hs .^ 2) # Stability¶

Are there “rough” functions for which these formulas estimate $$u'(x_i) = 0$$?

x = LinRange(-1, 1, 9)
f_rough(x) = cos(.1 + 4π*x)
fp_rough(x) = -4π*sin(.1 + 4π*x)

plot(x, f_rough, marker=:circle)
plot!(f_rough) fig = plot(fp_rough, xlims=(-1, 1))
for d in difflist
xx, yy = d(x, f_rough.(x))
plot!(fig, xx, yy, label=d, marker=:circle)
end
fig If we have a solution $$u(x)$$, then $$u(x) + f_{\text{rough}}(x)$$ is indistinguishable to our FD method.

# Second derivatives¶

We can compute a second derivative by applying first derivatives twice.

function diff2a(x, u)
xx, yy = diff1c(x, u)
diff1c(xx, yy)
end

function diff2b(x, u)
xx, yy = diff1l(x, u)
diff1r(xx, yy)
end

diff2list = [diff2a, diff2b]
n = 10
x = LinRange(-3, 3, n)
u = - cos.(x);

fig = plot(cos, xlims=(-3, 3))
for d2 in diff2list
xx, yy = d2(x, u)
plot!(fig, xx, yy, marker=:circle, label=d2)
end
fig # How fast do these approximations converge?¶

grids = 2 .^ (3:10)
hs = 1 ./ grids
function refinement_error2(f, f_xx, d2)
error = []
for n in grids
x = LinRange(-3, 3, n)
xx, yy = d2(x, f.(x))
push!(error, norm(yy - f_xx.(xx), Inf))
end
error
end

refinement_error2 (generic function with 1 method)

fig = plot(xscale=:log10, yscale=:log10)
for d2 in diff2list
error = refinement_error2(x -> -cos(x), cos, d2)
plot!(fig, hs, error, marker=:circle, label=d2)
end
plot!(fig, hs, hs .^ 2) # Differentiation matrices¶

All our diff* functions thus far have been linear in u, therefore they can be represented as matrices.

$\begin{split}\frac{u_{i+1} - u_i}{x_{i+1} - x_i} = \begin{bmatrix} -1/h & 1/h \end{bmatrix} \begin{bmatrix} u_i \\ u_{i+1} \end{bmatrix}\end{split}$

function diff1_mat(x)
n = length(x)
D = zeros(n, n)
h = x - x
D[1, 1:2] = [-1/h  1/h]
for i in 2:n-1
D[i, i-1:i+1] = [-1/2h  0  1/2h]
end
D[n, n-1:n] = [-1/h  1/h]
D
end

diff1_mat (generic function with 1 method)

x = LinRange(-3, 3, 10)
plot(x, diff1_mat(x) * sin.(x), marker=:circle)
plot!(cos) # How accurate is this derivative matrix?¶

fig = plot(xscale=:log10, yscale=:log10, legend=:topleft)
error = refinement_error(sin, cos, (x, u) -> (x, diff1_mat(x) * u))
plot!(fig, hs, error, marker=:circle)
plot!(fig, hs, hs, label="\$h\$")
plot!(fig, hs, hs .^ 2, label="\$h^2\$") # Can we study it as a matrix?¶

D = diff1_mat(x)
spy(D, marker=(:square, 10), c=:bwr) svdvals(D)

10-element Vector{Float64}:
2.268133218393964
2.2674392839412794
1.4265847744427302
1.368373796830966
1.2135254915624205
1.0228485194005286
0.8816778784387096
0.5437139466339259
0.46352549156242107
3.7873060138463766e-17


# Second derivative with Dirichlet boundary conditions¶

The left endpoint in our example boundary value problem 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{split} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ \\ & & A_{2:,:} & & \\ \\ \end{bmatrix} \begin{bmatrix} u_0 \\ \\ u_{2:} \\ \\ \end{bmatrix} = \begin{bmatrix} a \\ \\ f_{2:} \\ \\ \end{bmatrix} . \end{split}$

• This matrix is not symmetric even if $$A$$ is.

function laplacian_dirichlet(x)
n = length(x)
D = zeros(n, n)
h = x - x
D[1, 1] = 1
for i in 2:n-1
D[i, i-1:i+1] = (1/h^2) * [-1, 2, -1]
end
D[n, n] = 1
D
end

laplacian_dirichlet (generic function with 1 method)


# Laplacian as a matrix¶

L = laplacian_dirichlet(x)
spy(L, marker=(:square, 10), c=:bwr) svdvals(L)

10-element Vector{Float64}:
8.745632414048261
8.01252539624902
6.887027178840083
5.501492841921461
4.019212498940009
2.6202304303995274
1.510174898914501
0.9274512680423916
0.6309302491710911
0.23924866094931388


# Shortcomings of our previous methods¶

• Only second order accurate (at best)

• Worse than second order on non-uniform grids

• Worse than second order at Neumann boundaries

• Boundary conditions break symmetry

# Interpolation by Vandermonde matrices¶

We can compute a polynomial

$p(x) = c_0 + c_1 x + c_2 x^2 + \dotsb$

that assumes function values $$p(x_i) = u_i$$ by solving a linear system with the Vandermonde matrix.

$\begin{split} \underbrace{\begin{bmatrix} 1 & x_0 & x_0^2 & \dotsb \\ 1 & x_1 & x_1^2 & \dotsb \\ 1 & x_2 & x_2^2 & \dotsb \\ \vdots & & & \ddots \end{bmatrix}}_V \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} u_0 \\ u_1 \\ u_2 \\ \vdots \end{bmatrix} .\end{split}$
function vander(x, k=nothing)
if k === nothing
k = length(x)
end
V = ones(length(x), k)
for j = 2:k
V[:, j] = V[:, j-1] .* x
end
V
end

vander (generic function with 2 methods)

vander(LinRange(-1, 1, 5))

5×5 Matrix{Float64}:
1.0  -1.0  1.0   -1.0    1.0
1.0  -0.5  0.25  -0.125  0.0625
1.0   0.0  0.0    0.0    0.0
1.0   0.5  0.25   0.125  0.0625
1.0   1.0  1.0    1.0    1.0


# Fitting a polynomial¶

k = 4
x = LinRange(-2, 2, k)
u = sin.(x)
V = vander(x)
c = V \ u
scatter(x, u, label="\$u_i\$", legend=:topleft)
plot!(x -> (vander(x, k) * c)[1,1], label="\$p(x)\$")
plot!(sin, label=sin) # Differentiating¶

We’re given the coefficients $$c = V^{-1} u$$ of the polynomial

$p(x) = c_0 + c_1 x + c_2 x^2 + \dotsb.$
What is

(5)\begin{align} p(0) &= c_0 \\ p'(0) &= c_1 \\ p''(0) &= c_2 \cdot 2\\ p^{(k)}(0) &= c_k \cdot k! . \end{align}
function fdstencil1(source, target)
"first derivative stencil from source to target"
x = source .- target
V = vander(x)]
inv(V)[2, :]' # as a row vector
end
plot([z -> fdstencil1(x, z) * u, cos], xlims=(-3,3))
scatter!(x, 0*x, label="grid points")

syntax: "function" at In:1 expected "end", got "]"

Stacktrace:
 top-level scope
@ In:4
 eval
@ ./boot.jl:360 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


# Arbitrary order¶

function fdstencil(source, target, k)
"kth derivative stencil from source to target"
x = source .- target
V = vander(x)
rhs = zero(x)'
rhs[k+1] = factorial(k)
rhs / V
end
fdstencil(x, 0.5, 2)

1×4 adjoint(::Vector{Float64}) with eltype Float64:
0.0703125  0.351563  -0.914062  0.492188

plot([z -> fdstencil(x, z, 2) * u,
z -> -sin(z)]) ## We didn’t call inv(V); what’s up?¶

$p(0) = s_0^0 u_0 + s_1^0 u_1 + s_2^0 u_2 + \dotsb = e_0^T \underbrace{V^{-1} u}_c = \underbrace{e_0^T V^{-1}}_{s^0} u$

# Convergence order¶

hs = 2 .^ -LinRange(-4, 10, 10)
function diff_error(u, du, h; n, k, z=0)
x = LinRange(-h, h, n) .+ .5
fdstencil(x, z, k) * u.(x) - du.(z)
end
errors = [diff_error(sin, t -> -sin(t), h, n=5, k=2, z=.5+0.1*h)
for h in hs]
plot(hs, abs.(errors), marker=:circle)
plot!(h -> h^3, label="\$h^?\$", xscale=:log10, yscale=:log10) ## Observations¶

• When using $$n=3$$ points, we fit a polynomial of degree 2 and have error $$O(h^3)$$ for interpolation $$p(0)$$.

• Each derivative gives up one order of accuracy in general.

• Centered diff on uniform grids can have extra cancellation (superconvergence)

• The Vandermonde matrix is notoriously ill-conditioned with many points $$n$$. We recommend using a stable algorithm from Fornberg.

# High order discretization of the Laplacian¶

## The Poisson problem $$-u_{xx} = f$$ with boundary conditions¶

function poisson(x, spoints, forcing; left=(0, zero), right=(0, zero))
n = length(x)
L = zeros(n, n)
rhs = forcing.(x)
for i in 2:n-1
jleft = min(max(1, i-spoints÷2), n-spoints+1)
js = jleft : jleft + spoints - 1
L[i, js] = -fdstencil(x[js], x[i], 2)
end
L[1,1:spoints] = fdstencil(x[1:spoints], x, left)
L[n,n-spoints+1:n] = fdstencil(x[n-spoints+1:n], x[n], right)
rhs = left(x)
rhs[n] = right(x[n])
L, rhs
end

poisson (generic function with 1 method)

L, b = poisson(LinRange(-1, 1, 6), 3, zero, left=(1, zero))
L

6×6 Matrix{Float64}:
-3.75   5.0   -1.25   0.0           0.0           0.0
-6.25  12.5   -6.25   0.0           0.0           0.0
0.0   -6.25  12.5   -6.25          0.0           0.0
0.0    0.0   -6.25  12.5          -6.25          0.0
0.0    0.0    0.0   -6.25         12.5          -6.25
0.0    0.0    0.0    9.25186e-18  -2.31296e-16   1.0


# Method of manufactured solutions¶

## Problem: analytic solutions to PDEs are hard to find¶

Let’s choose a smooth function with rich derivatives,

$u(x) = \tanh(x) .$
Then
$u'(x) = \cosh^{-2}(x)$
and
$u''(x) = -2 \tanh(x) \cosh^{-2}(x) .$

• This works for nonlinear too.

x = LinRange(-2, 2, 10)
L, rhs = poisson(x, 3,
x -> 2 * tanh(x) / cosh(x)^2,
left=(0, tanh),
right=(1, x -> cosh(x)^-2))
u = L \ rhs
plot(x, u, marker=:circle, legend=:topleft)
plot!(tanh) # Convergence rate¶

ns = 2 .^ (2:10)
hs = 1 ./ ns
function poisson_error(n)
x = LinRange(-2, 2, n)
L, rhs = poisson(x, 3, x -> 2 * tanh(x) / cosh(x)^2,
left = (0, tanh),
right = (1, x -> cosh(x)^-2))
u = L \ rhs
norm(u - tanh.(x), Inf)
end

poisson_error (generic function with 1 method)

plot(hs, [poisson_error(n) for n in ns], marker=:circle)
plot!(h -> h^2, xscale=:log10, yscale=:log10) # Symmetry in boundary conditions: Dirichlet¶

We have implemented Dirichlet conditions by modifying the first row of the matrix,

$\begin{split} \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} . \end{split}$

• This matrix is not symmetric even if $$A$$ is.

• We can eliminate $$u_0$$ and create a reduced system for $$u_{1:}$$.

• Generalize: consider a $$2\times 2$$ block system

$\begin{split} \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} .\end{split}$

We can rearrange as

$A_{11} u_1 = f_1 - A_{10} f_0,$
which is symmetric if $$A_{11}$$ is.

• This is 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.

• Convenient alternative: write

$\begin{split} \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}, \end{split}$
which is symmetric and decouples the degrees of freedom associated with the boundary. This method applies cleanly to nonlinear problems.

• Optionally scale the identity by some scalar related to the norm of $$A_{11}$$.

# Symmetry in boundary conditions: Neumann¶

Consider FD discretization of the Neumann boundary condition

$\frac{du}{dx}(1) = b .$

1. Use a one-sided difference formula as in

$\frac{u_n - u_{n-1}}{h} = b .$

• an extra discretization choice

• may reduce order of accuracy compared to interior discretization, lose symmetry.

1. 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,

$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.

# Fourier analysis of stencils¶

Consider the plane waves $$\phi(x, \theta) = e^{i\theta x}$$.

Sample $$\phi$$ on a discrete grid $$x = \mathbb Z$$ and apply the stencil

(6)\begin{align} S \phi(x, \theta) &= s_{-1} \phi(x-1, \theta) + s_{0} \phi(x, \theta) + s_1 \phi(x+1,\theta) \\ &= \Big( s_{-1} e^{-i\theta} + s_0 + s_{1} e^{i\theta} \Big) \phi(x, \theta) \end{align}

With $$S = \begin{bmatrix} -1 & 2 & -1 \end{bmatrix}$$, we get

$S \phi(x, \theta) = \underbrace{(2 - 2 \cos\theta)}_{\hat S(\theta)} \phi(x, \theta)$
We call $$\hat S(\theta)$$ the symbol of the operator. What is the symbol of the continuous second derivative?

plot(theta -> 2 - 2*cos(theta), xlims=(-pi, pi)) 