.ipynb .pdf

repository open issue

Binder Live Code

Contents

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[2] - x[1]
    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[2] - x[1]
    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[23]:1 expected "end", got "]"

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

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

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

previous

2021-08-27 Finite Difference 2

next

2021-09-01 Boundaries and Fourier methods