2022-09-02 Chebyshev collocation#

Last time#

Solutions and matrix properties
Method of Manufactured Solutions
Techniques for boundary conditions
Fourier analysis of stencils

Today#

Conditioning of Vandermonde matrices
Chebyshev polynomials
Solving Poisson using Chebyshev polynomials

using Plots
default(linewidth=3)
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_fd(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_fd (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,

\[ 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, 15)
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)

UndefVarError: poisson not defined

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

Convergence rate#

ns = 2 .^ (3:10)
hs = 1 ./ ns
function poisson_fd_error(n; spoints=3)
    x = LinRange(-2, 2, n)
    L, rhs = poisson_fd(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_fd_error (generic function with 1 method)

plot(hs, [poisson_fd_error(n, spoints=5) for n in ns], marker=:circle)
plot!(h -> h^4, label="\$h^2\$", xscale=:log10, yscale=:log10)

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

theta = LinRange(-pi, pi, 10)
symbol([-.5 0 .5], theta)
#symbol([1 -2 1], theta)

10-element Vector{ComplexF64}:
0 - 1.2246467991473532e-16im
0 - 0.6427876096865395im
0 - 0.9848077530122081im
0 - 0.8660254037844388im
0 - 0.34202014332566877im
0 + 0.34202014332566877im
0 + 0.8660254037844385im
0 + 0.9848077530122081im
0 + 0.6427876096865395im
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)

../_images/2022-09-02-chebyshev_11_0.svg

Stencils of high-order operators#

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

../_images/2022-09-02-chebyshev_13_0.svg

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

../_images/2022-09-02-chebyshev_14_0.svg

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

Visualizing matrix transformations#

function peanut()
    theta = LinRange(0, 2*pi, 50)
    r = 1 .+ .4*sin.(3*theta) + .6*sin.(2*theta)
    x = r .* cos.(theta)
    y = r .* sin.(theta)
    x, y
end
x, y = peanut()
scatter(x, y, aspect_ratio=:equal)

../_images/2022-09-02-chebyshev_17_0.svg

Group these points into a \(2\times n\) matrix \(X\). Note that multiplication by any matrix \(A\) is applied to each column separately, i.e.,

\[ A \underbrace{\Bigg[ \mathbf x_1 \Bigg| \mathbf x_2 \Bigg| \dotsb \Bigg]}_X = \Bigg[ A \mathbf x_1 \Bigg| A \mathbf x_2 \Bigg| \dotsb \Bigg] \]

X = [x y]'
size(X)

(2, 50)

Visualizing the Singular Value Decomposition#

function tplot(A)
    x, y = peanut()
    X = [x y]'
    Y = A * X
    scatter(X[1,:], X[2,:], label="input", aspect_ratio=:equal)
    scatter!(Y[1,:], Y[2,:], label="output")
end

A = randn(2, 2)
display(svdvals(A))
tplot(A)

2-element Vector{Float64}:
 1.1456971880440263
 0.06488883506970833

../_images/2022-09-02-chebyshev_21_1.svg

U, S, V = svd(A)
tplot(U)

../_images/2022-09-02-chebyshev_22_0.svg

Condition number of interpolation#

n = 20
source = LinRange(-1, 1, n)
target = LinRange(-1, 1, 100)
P = vander(source)
A = vander(target, n) / P
svdvals(A)

20-element Vector{Float64}:
821782616643
1886116478436
201394442903565
2003668113747805
326389915733187
28865315753632
2826626548042968
282658231256019
282657730773852
2826577307583236
2826577297323856
2826577141906412
2825938632428557
282352200333904
2484549739833892
2449230244092657
9014541632841195
7794033224148031
0093167855491605
0064482710718898

runge(x) = 1 / (1 + 10*x^2)
y = runge.(source)
scatter(source, y)
plot!([runge,
        x -> (vander([x], n) / P * y)[1]])
# plot!(target, A * y)

../_images/2022-09-02-chebyshev_25_0.svg

The bad singular vectors#

\[ A = U \Sigma V^T \]

n = 11
source = LinRange(-1, 1, n)
target = LinRange(-1, 1, 100)
P = vander(source)
A = vander(target, n) / P
U, S, V = svd(A)
S

11-element Vector{Float64}:
77706851205524
13125653860779
200622823579756
1551891145187843
1464270560667247
146390542844354
1459494752514874
900733962570026
700731526422623
2170432449283064
1628054393968887

scatter(source, V[:,1])
plot!(target, U[:,1])

../_images/2022-09-02-chebyshev_28_0.svg

Choosing the source points#

CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))

CosRange (generic function with 1 method)

n = 15
source = CosRange(-1, 1, n)
target = LinRange(-1, 1, 100)
P = vander(source)
A = vander(target, n) / P
svdvals(A)

15-element Vector{Float64}:
322445257960405
295148962438826
239022400818026
1811090404013336
0703313983108265
9872800167943647
812157783053922
705974980411089
4581412991002862
320177305779112
017320105389942
7903748104690793
5900555135944288
0373676954023048
031069620176098

runge(x) = 1 / (1 + 10*x^2)
y = runge.(source)
scatter(source, y)
plot!([runge,
        x -> (vander([x], n) / P * y)[1]])

../_images/2022-09-02-chebyshev_32_0.svg

Condition number of interpolation on `CosRange` points#

interpolate(source, target=LinRange(-1, 1, 100)) = vander(target, length(source)) / vander(source)

plot_cond(mat, points) = plot!([
        cond(mat(points(-1, 1, n)))
        for n in 2:30], label="$mat/$points", marker=:auto, yscale=:log10)

plot(title="Interpolation condition numbers")
plot_cond(interpolate, LinRange)
plot_cond(interpolate, CosRange)

../_images/2022-09-02-chebyshev_34_0.svg

Vandermonde conditioning, stable algorithms#

It is possible for interpolation to be well-conditioned, but construct it from ill-conditioned pieces.

plot(title="Vandermonde condition numbers")
plot_cond(vander, LinRange)
plot_cond(vander, CosRange)

../_images/2022-09-02-chebyshev_36_0.svg

Chebyshev polynomials#

function vander_chebyshev(x, n=nothing)
    if isnothing(n)
        n = length(x) # Square by default
    end
    m = length(x)
    T = ones(m, n)
    if n > 1
        T[:, 2] = x
    end
    for k in 3:n
        T[:, k] = 2 * x .* T[:,k-1] - T[:, k-2]
    end
    T
end

vander_chebyshev (generic function with 2 methods)

plot(title="Vandermonde condition numbers")
plot_cond(vander, LinRange)
plot_cond(vander, CosRange)
plot_cond(vander_chebyshev, LinRange)
plot_cond(vander_chebyshev, CosRange)

../_images/2022-09-02-chebyshev_39_0.svg

Derivation of Chebyshev polynomials (supplemental)#

Chebyshev polynomials are defined for non-negative integar \(n\) as

\[ T_n(x) = \cos (n \arccos(x)) .\]

This turns out to be a polynomial, but it’s not obvious why. Recall

\[ \cos(a + b) = \cos a \cos b - \sin a \sin b .\]

Let \(y = \arccos x\) and check

\[\begin{split} \begin{split} T_{n+1}(x) &= \cos \big( (n+1) y \big) = \cos ny \cos y - \sin ny \sin y \\ T_{n-1}(x) &= \cos \big( (n-1) y \big) = \cos ny \cos y + \sin ny \sin y \end{split}\end{split}\]

Adding these together produces

\[ T_{n+1}(x) + T_{n-1}(x) = 2\cos ny \cos y = 2 x \cos ny = 2 x T_n(x) \]

which yields a convenient recurrence:

\[\begin{split}\begin{split} T_0(x) &= 1 \\ T_1(x) &= x \\ T_{n+1}(x) &= 2 x T_n(x) - T_{n-1}(x) \end{split}\end{split}\]

Lagrange interpolating polynomials#

Find the unique polynomial that is 0 at all but one source point.

source = LinRange(-1, 1, 11)
target = LinRange(-1, 1, 200)
scatter(source, one.(source), title="LinRange")
plot!(target, interpolate(source, target))

../_images/2022-09-02-chebyshev_42_0.svg

source = CosRange(-1, 1, 13)
target = LinRange(-1, 1, 200)
scatter(source, one.(source), title="CosRange")
plot!(target, interpolate(source, target))

../_images/2022-09-02-chebyshev_43_0.svg

Differentiating via Chebyshev#

We can differentiate Chebyshev polynomials using the recurrence

\[ \frac{T_k'(x)}{k} = 2 T_{k-1}(x) + \frac{T_{k-2}'(x)}{k-2} \]

which we can differentiate to evaluate higher derivatives.

function chebdiff(x, n=nothing)
    T = vander_chebyshev(x, n)
    m, n = size(T)
    dT = zero(T)
    dT[:,2:3] = [one.(x) 4*x]
    for j in 3:n-1
        dT[:,j+1] = j * (2 * T[:,j] + dT[:,j-1] / (j-2))
    end
    ddT = zero(T)
    ddT[:,3] .= 4
    for j in 3:n-1
        ddT[:,j+1] = j * (2 * dT[:,j] + ddT[:,j-1] / (j-2))
    end
    T, dT, ddT
end

chebdiff (generic function with 2 methods)

x = CosRange(-1, 1, 20)
T, dT, ddT = chebdiff(x)
c = T \ cos.(3x)
scatter(x, dT * c)
plot!(x -> -3sin(3x))

../_images/2022-09-02-chebyshev_46_0.svg

Solving a BVP on Chebyshev nodes#

function poisson_cheb(n, rhsfunc, leftbc=(0, zero), rightbc=(0, zero))
    x = CosRange(-1, 1, n)
    T, dT, ddT = chebdiff(x)
    L = -ddT
    rhs = rhsfunc.(x)
    for (index, deriv, func) in
            [(1, leftbc...), (n, rightbc...)]
        L[index,:] = (T, dT)[deriv+1][index,:]
        rhs[index] = func(x[index])
    end
    x, L / T, rhs
end

poisson_cheb (generic function with 3 methods)

manufactured(x) = tanh(2x)
d_manufactured(x) = 2*cosh(2x)^-2
mdd_manufactured(x) = 8 * tanh(2x) / cosh(2x)^2
x, A, rhs = poisson_cheb(12, mdd_manufactured,
    (0, manufactured), (1, d_manufactured))
plot(x, A \ rhs, marker=:auto)
plot!(manufactured, legend=:bottomright)

../_images/2022-09-02-chebyshev_49_0.svg

“spectral” (exponential) convergence#

function poisson_error(n)
    x, A, rhs = poisson_cheb(n, mdd_manufactured, (0, manufactured), (1, d_manufactured))
    u = A \ rhs
    norm(u - manufactured.(x), Inf)
end

ns = 3:40
plot(ns, abs.(poisson_error.(ns)), marker=:circle, yscale=:log10)
ps = [1 2 3]
plot!([n -> n^-p for p in ps], label=map(p -> "\$n^{-$p}\$", ps))

../_images/2022-09-02-chebyshev_51_0.svg

plot(ns, [poisson_fd_error(n, spoints=3) for n in ns], 
    marker=:circle, yscale=:log10)
plot!([n -> n^-p for p in ps], label=map(p -> "\$n^{-$p}\$", ps))

../_images/2022-09-02-chebyshev_52_0.svg

Numerical Solution of Partial Differential Equations

2022-09-02 Chebyshev collocation

Contents

2022-09-02 Chebyshev collocation#

Last time#

Today#

Method of manufactured solutions#

Problem: analytic solutions to PDEs are hard to find#

Convergence rate#

Numerically computing symbols#

Stencils of high-order operators#

Outlook on Fourier methods#

Visualizing matrix transformations#

Visualizing the Singular Value Decomposition#

Condition number of interpolation#

The bad singular vectors#

Choosing the source points#

Condition number of interpolation on `CosRange` points#

Vandermonde conditioning, stable algorithms#

Chebyshev polynomials#

Derivation of Chebyshev polynomials (supplemental)#

Lagrange interpolating polynomials#

Differentiating via Chebyshev#

Solving a BVP on Chebyshev nodes#

“spectral” (exponential) convergence#

Numerical Solution of Partial Differential Equations

2022-09-02 Chebyshev collocation

Contents

2022-09-02 Chebyshev collocation#

Last time#

Today#

Method of manufactured solutions#

Problem: analytic solutions to PDEs are hard to find#

Convergence rate#

Numerically computing symbols#

Stencils of high-order operators#

Outlook on Fourier methods#

Visualizing matrix transformations#

Visualizing the Singular Value Decomposition#

Condition number of interpolation#

The bad singular vectors#

Choosing the source points#

Condition number of interpolation on CosRange points#

Vandermonde conditioning, stable algorithms#

Chebyshev polynomials#

Derivation of Chebyshev polynomials (supplemental)#

Lagrange interpolating polynomials#

Differentiating via Chebyshev#

Solving a BVP on Chebyshev nodes#

“spectral” (exponential) convergence#

Condition number of interpolation on `CosRange` points#