# 2021-09-01 Boundaries and Fourier methods¶

## Last time¶

• High order on arbitrary grids

• Method of manufactured solutions

• More satisfying approach to boundary conditions

## Today¶

• More on boundary conditions

• Fourier methods for analyzing finite difference methods

• Discussion

using Plots
using LinearAlgebra


# 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)

cond(vander(LinRange(-1, 1, 5) .+ 1))
#cond(vander([-1,-.99,-.98,-.97, 1]))

743.2516243247674


# Arbitrary order derivatives¶

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
x = LinRange(-2, 2, 5)
u = sin.(x)
fdstencil(x, 0.5, 2)

1×5 adjoint(::Vector{Float64}) with eltype Float64:
-0.208333  1.33333  -1.75  0.333333  0.291667

plot([z -> fdstencil(x, z, 2) * u,
z -> -sin(z)]) # 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, 5,
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 .^ (3:10)
hs = 1 ./ ns
function poisson_error(n, spoints)
x = LinRange(-2, 2, n)
L, rhs = poisson(x, spoints, 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, 5) for n in ns], marker=:circle)
plot!(hs, [poisson_error(n, 3) for n in ns], marker=:square)
plot!([h -> h^2, h -> h^4], xscale=:log10, yscale=:log10, legend=:topleft) # Symmetry in boundary conditions: Dirichlet¶

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

(7)\begin{align} \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{align}
• 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{align*} \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{align*}

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.