# 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


```julia
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
```

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


```julia
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)


```julia
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)
```




    

    




```julia
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)


```julia
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)
```




    

    




```julia
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


```julia
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)




```julia
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 solution
  * manufactured 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](https://doi.org/10.1115/1.1436090) (7 pages)

# Symmetry in boundary conditions: Dirichlet

We have implemented Dirichlet conditions 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} . $$

* 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{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 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{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 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.
2. 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
\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?


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




    

    



# Numerically computing symbols


```julia
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




```julia
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


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




    

    




```julia
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

76a14134ab51716d0e28f4b9762ccd663681f3c5
ipynb
Jupyter Notebook
slides/2021-09-03-fourier.ipynb
cu-numpde/numpde
e5e1a465a622eba56900004f9a503412407cdccf
MIT
2022-01-01T20:54:51.000Z
2022-01-01T20:54:51.000Z
slides/2021-09-03-fourier.ipynb
amta3208/fall21
e5e1a465a622eba56900004f9a503412407cdccf
MIT
slides/2021-09-03-fourier.ipynb
amta3208/fall21
e5e1a465a622eba56900004f9a503412407cdccf
MIT
2022-01-01T20:54:46.000Z
2022-01-01T20:54:46.000Z
Qwen/Qwen-72B
1. YES
2. YES