# 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


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

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


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




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




    743.2516243247674



# Arbitrary order derivatives


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




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


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




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


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


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




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




    

    



# Symmetry in boundary conditions: Dirichlet

We have implemented Dirichlet conditions by modifying the first row of the matrix,
\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{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.

a01db38d0d98e13bb6c51b69ef841ea9d54edfad
ipynb
Jupyter Notebook
slides/2021-09-01-convergence-and-boundaries.ipynb
cu-numpde/numpde
e5e1a465a622eba56900004f9a503412407cdccf
MIT
2022-01-01T20:54:51.000Z
2022-01-01T20:54:51.000Z
slides/2021-09-01-convergence-and-boundaries.ipynb
amta3208/fall21
e5e1a465a622eba56900004f9a503412407cdccf
MIT
slides/2021-09-01-convergence-and-boundaries.ipynb
amta3208/fall21
e5e1a465a622eba56900004f9a503412407cdccf
MIT
2022-01-01T20:54:46.000Z
2022-01-01T20:54:46.000Z
Qwen/Qwen-72B
1. YES
2. YES