2023-09-01 FD Solutions

2023-09-01 FD Solutions#

Last time#

Discussion and scoping
Evaluating derivatives
Taylor series and truncation error

Today#

Discuss activity
Measuring errors
Stable discretizations
Solutions and matrix properties

using Plots
default(linewidth=3)
using LinearAlgebra

Finite Difference/collocation approach to solve \(u\):#

(5)#\[\begin{gather} -\frac{d^2 u}{dx^2} = f(x) \quad x \in \Omega = (-1,1) \\ u(-1) = a \quad \frac{du}{dx}(1) = b . \end{gather}\]

Represent the function \(u(x)\) by its values \(u_i = u(x_i)\) at a discrete set of points
\[ -1 = x_1 < x_2 < \dotsb < x_n = 1 . \]
- The FD framework does not uniquely specify the solution values at other points
Compute derivatives at \(x_i\) via differencing formulas involving a finite number of neighbor points (independent of the total number of points \(n\)).
FD methods ask for the differential equation to be satisfied pointwise at each \(x_i\) in the interior of the domain.
Evaluate the forcing term \(f\) pointwise at \(x_i\).
Approximate derivatives at discrete boundary points (\(x_n = 1\) above), typically using one-sided differencing formulas.

Activity on stability and accuracy computing derivatives#

A note on stable computation#

x = big(1e-15)
@show big(1+x)
@show log(1 + x)

big(1 + x) = 1.000000000000001000000000000000077705399876661079238307185601195015145492561714
log(1 + x) = 9.999999999999995777053998766613348662406422732737965069319591605458493216045864e-16

9.999999999999995777053998766613348662406422732737965069319591605458493216045864e-16

plot([h -> log(1+h), log1p], xlim=(-1e-15, 1e-15))

What happens as we zoom in?

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 = 40 
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
default(legendfontsize=12)

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), 2)/sqrt(n))
    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 hs .^ 2], label=["h" "\$h^2\$"])

What happens if we use a 1-norm or 2-norm?

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.

Consistency + Stability = Convergence#

[Lax Equivalence Theorem]#

Consistency#

When we apply the differential operator to the exact solution, we get a small residual.
The residual converges under grid refinement.
Hopefully fast as \(h \to 0\)

Stability#

There do not exist “bad” functions that also satisfy the equation.
The solution does not “blow up” for time-dependent problems.
Definition here is intentionally vague, and there are more subtle requirements for problems like incompressible flow.

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 = 20
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
        )) # which norm?
    end
    error
end

refinement_error2 (generic function with 1 method)

fig = plot(xlabel="h", xscale=:log10, ylabel="Error", 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, label="\$h^2\$") 

Both methods are second order accurate.
The diff2b method is more accurate than diff2a (by a factor of 4)
The diff2a method can’t compute derivatives at points adjacent the boundary.
We don’t know yet whether either is stable

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
x = LinRange(-1, 1, 5)
diff1_mat(x)

5×5 Matrix{Float64}:
 -2.0   2.0   0.0   0.0  0.0
 -1.0   0.0   1.0   0.0  0.0
  0.0  -1.0   0.0   1.0  0.0
  0.0   0.0  -1.0   0.0  1.0
  0.0   0.0   0.0  -2.0  2.0

n = 12
x = LinRange(-3, 3, n)
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?#

function my_spy(A)
    cmax = norm(vec(A), Inf)
    s = max(1, ceil(120 / size(A, 1)))
    spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

D = diff1_mat(x)
my_spy(D)

svdvals(D)

12-element Vector{Float64}:
7718000318790152
7716788692501084
7708640148632901
732003875931066
5877132402714698
4519108222799693
2963624321753364
0397961559479483
9166666666666657
5413728320040285
4745015826879543
7546012240164312e-16

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_1 = u(x_1 = -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_1 \\ \\ 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)
my_spy(L)

cond(L)

56.42029373680045

Solutions#

\[\begin{split} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ \\ & & A_{2:n-1,:} & & \\ \\ \\ 0&0 &0 &0 &0 & 1 \end{bmatrix} \begin{bmatrix} u_1 \\ \\ u_{2:} \\ \\ \\ u_n \end{bmatrix} = \begin{bmatrix} a \\ \\ f_{2:} \\ \\ \\ b \end{bmatrix} . \end{split}\]

L = laplacian_dirichlet(x)
f = one.(x)
f[1] = 0
f[end] = 0;
plot(x, f)

u = L \ f
plot(x, u)

Discrete “Green’s functions”#

plot(x, inv(L)[:, 2])

Ln = copy(L)
Ln[n, n-1] = -1
plot(x, inv(Ln)[:, 4])

Discrete eigenfunctions#

x = LinRange(-1, 1, 10)
L = laplacian_dirichlet(x)
Lambda, V = eigen(L)
plot(Lambda, marker=:circle)

plot(x, V[:, 1:4])

2023-09-01 FD Solutions

Contents

2023-09-01 FD Solutions#

Last time#

Today#

Finite Difference/collocation approach to solve \(u\):#

Activity on stability and accuracy computing derivatives#

A note on stable computation#

A few methods on grids#

Measuring error on grids#

Stability#

Consistency + Stability = Convergence#

[Lax Equivalence Theorem]#

Consistency#

Stability#

Second derivatives#

How fast do these approximations converge?#

Differentiation matrices#

How accurate is this derivative matrix?#

Can we study it as a matrix?#

Second derivative with Dirichlet boundary conditions#

Laplacian as a matrix#

Solutions#

Discrete “Green’s functions”#

Discrete eigenfunctions#

Outlook on our method#

Pros#

Cons#