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


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

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


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

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


# 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}:
4141.821782616643
756.1886116478436
11.201394442903565
5.2003668113747805
2.326389915733187
2.28865315753632
2.2826626548042968
2.282658231256019
2.282657730773852
2.2826577307583236
2.2826577297323856
2.2826577141906412
2.2825938632428557
2.282352200333904
2.2484549739833892
2.2449230244092657
1.9014541632841195
1.7794033224148031
1.0093167855491605
1.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)


$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}:
34.77706851205524
11.13125653860779
3.200622823579756
3.1551891145187843
3.1464270560667247
3.146390542844354
3.1459494752514874
2.900733962570026
2.700731526422623
1.2170432449283064
1.1628054393968887

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


# 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}:
3.322445257960405
3.295148962438826
3.239022400818026
3.1811090404013336
3.0703313983108265
2.9872800167943647
2.812157783053922
2.705974980411089
2.4581412991002862
2.320177305779112
2.017320105389942
1.7903748104690793
1.5900555135944288
1.0373676954023048
1.031069620176098

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


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


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


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


# 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}$
$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))

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


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


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


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