# 2021-09-03 Convergence and Fourier analysis¶

## Last time¶

Method of manufactured solutions

More on boundary conditions and symmetry

## Today¶

Investigate convergence with manufactured solutions

Fourier methods for analyzing finite difference methods

Discussion and introduce/start boundary conditions activity

```
using Plots
using LinearAlgebra
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
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
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[1], left[1])
L[n,n-spoints+1:n] = fdstencil(x[n-spoints+1:n], x[n], right[1])
rhs[1] = left[2](x[1])
rhs[n] = right[2](x[n])
L, rhs
end
```

```
poisson (generic function with 1 method)
```

# Method of manufactured solutions¶

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

Let’s choose a smooth function with rich derivatives,

This works for nonlinear too.

```
x = LinRange(-2, 2, 5)
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)
```

# Compare MMS error with different number of stencil points (Neumann)¶

```
x = LinRange(-2, 2, 5)
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)
```

```
x = LinRange(-2, 2, 5)
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)
```

# Compare MMS error with different number of stencil points (Dirichlet)¶

```
x = LinRange(-2, 2, 5)
L, rhs = poisson(x, 3,
x -> 2 * tanh(x) / cosh(x)^2,
left=(0, tanh),
right=(0, tanh))
u = L \ rhs
plot(x, u, marker=:circle, legend=:topleft)
plot!(tanh)
```

```
x = LinRange(-2, 2, 5)
L, rhs = poisson(x, 5,
x -> 2 * tanh(x) / cosh(x)^2,
left=(0, tanh),
right=(0, tanh))
u = L \ rhs
plot(x, u, marker=:circle, legend=:topleft)
plot!(tanh)
```

# Convergence rate¶

```
ns = [5, 6, 8, 10, 15, 30, 100]
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)
```

## Discuss¶

asymptotic vs preasmptotic regime

errors on coarse grids (here especially) may be dominated by the boundary conditions

leading error coeffficient in Taylor expansion is bigger near endpoints

model can be considered

**verified**when it converges at the design order of accuracy for one “rich” manufactured solutionmanufactured solution needs to exercise every aspect of the model

all types/cases of boundary conditions

a reliable tool only for smooth problems

for more: Roache, Code Verification by the Method of Manufactured Solutions (7 pages)

# Symmetry in boundary conditions: Dirichlet¶

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

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

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

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.

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

which simplifies to

# 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

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

*symbol*of the operator. What is the symbol of the continuous second derivative?

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

# Numerically computing symbols¶

```
function symbol(S, theta)
if length(S) % 2 != 1
error("Length of stencil must be odd")
end
w = length(S) ÷ 2
phi = exp.(1im * (-w:w) * theta')
vec(S * phi) # not! (S * phi)'
end
symbol([-.5 0 .5], LinRange(-pi, pi, 10))
```

```
10-element Vector{ComplexF64}:
0.0 - 1.2246467991473532e-16im
0.0 - 0.6427876096865395im
0.0 - 0.9848077530122081im
0.0 - 0.8660254037844388im
0.0 - 0.34202014332566877im
0.0 + 0.34202014332566877im
0.0 + 0.8660254037844385im
0.0 + 0.9848077530122081im
0.0 + 0.6427876096865395im
0.0 + 1.2246467991473532e-16im
```

```
function plot_symbol(S, deriv, n_theta=30)
theta = LinRange(-pi, pi, n_theta)
sym = symbol(S, theta)
rsym = real.((-1im)^deriv * sym)
fig = plot(theta, rsym, marker=:circle, label="stencil")
plot!(fig, th -> th^deriv, label="exact")
fig
end
plot_symbol([-.5 0 .5], 1)
#plot_symbol([1 -2 1], 2)
```

# Stencils of high-order operators¶

```
x = -5:5
plot_symbol(fdstencil(x, 0, 1), 1)
```

```
x = -5:5
plot_symbol(fdstencil(x, 0, 2), 2)
```

# Outlook on Fourier methods¶

the Fourier modes \(e^{i\theta x}\) and their multi-dimensional extensions are eigenvectors of

**all**stencil-type operations“high frequencies” \([-\pi, \pi) \setminus [-\pi/2, \pi/2)\) are generally poorly resolved so we need to use a grid fine enough that important features are at low frequencies \([-\pi/2, \pi/2)\)

same technique can be used to study the inverse (and approximations thereof), as with multigrid and multilevel domain decomposition methods (later in the course)

these methods can also be framed within the theory of (block) Toeplitz matrices