Chapter 13 - Fundamentals of Numerical Computation

Functions¶

Create a tensor-product grid

tensorgrid.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
"""
    tensorgrid(x, y)

Create a tensor grid for a rectangle from its 1d projections `x` and `y`.
Returns `unvec` to reshape a vector into a 2d array, `mtx` to evaluate a
function on the grid, `X`, `Y` to give the grid coordinates, and boolean array
`is_boundary` to identify the boundary points.
"""
function tensorgrid(x, y)
    m, n = length(x) - 1, length(y) - 1
    unvec(u) = reshape(u, m+1, n+1)
    mtx(h) = [h(x, y) for x in x, y in y]
    X = mtx((x, y) -> x)
    Y = mtx((x, y) -> y)
    is_boundary = trues(m+1, n+1)
    is_boundary[2:m, 2:n] .= false
    return mtx, X, Y, unvec, is_boundary
end

Solution of Poisson’s equation by finite differences

poissonfd.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
"""
    poissonfd(f, g, m, xspan, n, yspan)

Solve Poisson's equation on a rectangle by finite differences.
Function `f` is the forcing function and function `g` gives the
Dirichlet boundary condition. The rectangle is the tensor product of
intervals `xspan` and `yspan`,  and the discretization uses `m`+1
and `n`+1 points in the two coordinates.

Returns vectors defining the grid and a matrix of grid solution values.
"""
function poissonfd(f, g, m, xspan, n, yspan)
    # Discretize the domain.
    x, Dx, Dxx = FNC.diffmat2(m, xspan)
    y, Dy, Dyy = FNC.diffmat2(n, yspan)
    mtx, X, Y, unvec, is_boundary = tensorgrid(x, y)
    N = (m+1) * (n+1)   # total number of unknowns

    # Form the collocated PDE as a linear system.
    A = kron(I(n+1), sparse(Dxx)) + kron(sparse(Dyy), I(m+1))
    b = vec(mtx(f))

    # Apply Dirichlet condition.
    scale = maximum(abs, A[n+2, :])
    idx = vec(is_boundary)
    A[idx, :] = scale * I(N)[idx, :]        # Dirichet assignment
    b[idx] = scale * g.(X[idx], Y[idx])    # assigned values

    # Solve the linear system and reshape the output.
    u = A \ b
    return x, y, unvec(u)
end

Solution of elliptic PDE by Chebyshev collocation

elliptic.jl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
"""
    elliptic(ϕ, g, m, xspan, n, yspan)

Solve the elliptic PDE
    `ϕ`(x, y, u, u_x, u_xx, u_y, u_yy) = 0
on the rectangle `xspan`x`yspan`, subject to `g`(x,y)=0 on the boundary. Uses
`m`+1 points in x by `n`+1 points in y in a Chebyshev discretization. Returns
vectors defining the grid and a matrix of grid solution values.
"""
function elliptic(ϕ, g, m, xspan, n, yspan)
    # Discretize the domain.
    x, Dx, Dxx = diffcheb(m, xspan)
    y, Dy, Dyy = diffcheb(n, yspan)
    mtx, X, Y, unvec, is_boundary = tensorgrid(x, y)
    N = (m+1) * (n+1)   # total number of unknowns

    # Identify boundary locations and evaluate the boundary condition.
    idx = vec(is_boundary)
    gb = g.(X[idx], Y[idx])

    # Evaluate the PDE+BC residual.
    function residual(u)
        U = unvec(u)
        R = ϕ(X, Y, U, Dx * U, Dxx * U, U * Dy', U * Dyy')    # PDE
        @. R[idx] = u[idx] - gb                               # boundary residual
        return vec(R)
    end

    # Solve the equation.
    u = levenberg(residual, vec(zeros(size(X))))[end]
    U = unvec(u)

    return function (ξ, η)
        v = [chebinterp(x, u, ξ) for u in eachcol(U)]
        return chebinterp(y, v, η)
    end
end

Examples¶

13.1 Tensor-product discretizations¶

Example 13.1.2

Here is the grid from Example 13.1.1.

m = 4
x = range(0, 2, m+1)
n = 2
y = range(1, 3, n+1);

For a given $f(x,y)$ we can find $\operatorname{mtx}(f)$ by using a comprehension syntax.

f = (x, y) -> cos(π * x * y - y)
F = [f(x, y) for x in x, y in y]

5×3 Matrix{Float64}:
  0.540302  -0.416147  -0.989992
  0.841471   0.416147  -0.14112
 -0.540302  -0.416147   0.989992
 -0.841471   0.416147   0.14112
  0.540302  -0.416147  -0.989992

We can make a nice plot of the function by first choosing a much finer grid. However, the contour and surface plotting functions expect the transpose of mtx( $f$ ).

using Plots
m, n = 80, 60
x = range(0, 2, m+1);
y = range(1, 3, n+1);
F = [f(x, y) for x in x, y in y]
contour(x, y, F';
    levels=21,  aspect_ratio=1,
    color=:redsblues,  clims=(-1, 1),
    xlabel="x",  ylabel="y" )

Example 13.1.3

For a function given in polar form, such as $f(r,\theta)=1-r^4$ , construction of a function over the unit disk is straightforward using a grid in $(r,\theta)$ space.

r = range(0, 1, 41)
θ = range(0, 2π, 81)
F = [1 - r^4 for r in r, θ in θ]
plot(r, θ, F';
    legend=:none, 
    color=:viridis,  fill=true,
    xlabel="r",  ylabel="θ", 
    title="A polar function")

Of course, we are used to seeing such plots over the $(x,y)$ plane, not the $(r,\theta)$ plane.

In such functions the values along the line $r=0$ must be identical, and the values on the line $\theta=0$ should be identical to those on $\theta=2\pi$ . Otherwise the interpretation of the domain as the unit disk is nonsensical. If the function is defined in terms of $x$ and $y$ , then those can be defined in terms of $r$ and θ using (13.1.6).

Example 13.1.4

We define a function and, for reference, its two exact partial derivatives.

u = (x, y) -> sin(π * x * y - y);
∂u_∂x = (x, y) -> π * y * cos(πx * y - y);
∂u_∂y = (x, y) -> (π * x - 1) * cos(π * x * y - y);

We use an equispaced grid and second-order finite differences as implemented by diffmat2.

m = 80;
x, Dx, _ = FNC.diffmat2(m, [0, 2]);
n = 60;
y, Dy, _ = FNC.diffmat2(n, [1, 3]);
mtx = (f, x, y) -> [f(x, y) for x in x, y in y]
U = mtx(u, x, y)
∂xU = Dx * U
∂yU = U * Dy';

Now we compare the exact $\frac{\partial u}{\partial y}$ with its finite-difference approximation.

M = maximum(abs, ∂yU)    # find the range of the result
plot(layout=(1, 2), 
    aspect_ratio=1,   clims=(-M, M), 
    xlabel="x", ylabel="y")
contour!(x, y, mtx(∂u_∂y, x, y)';
    levels=15,  subplot=1,
    color=:redsblues,
    title="∂u/∂y")
contour!(x, y, ∂yU';
    levels=15,  subplot=2,
    color=:redsblues, 
    title="approximation")

To the eye there is little difference to be seen, though the results have no more than a few correct digits at these discretization sizes:

exact = mtx(∂u_∂y, x, y)
# Relative difference in Frobenius norm:
norm(exact - ∂yU) / norm(exact)

0.0035544848411698123

13.2 Two-dimensional diffusion and advection¶

Example 13.2.1

m = 2;
n = 3;
V = rand(1:9, m, n);
v = vec(V)

6-element Vector{Int64}:
 8
 4
 3
 1
 2
 9

The unvec operation is the inverse of vec.

unvec = z -> reshape(z, m, n)
unvec(v)

2×3 Matrix{Int64}:
 8  3  2
 4  1  9

Example 13.2.2

m, n = (60, 25)
x, Dx, Dxx = FNC.diffper(m, [-1, 1])
y, Dy, Dyy = FNC.diffper(n, [-1, 1])
mtx = f -> [f(x, y) for x in x, y in y]
unvec = z -> reshape(z, m, n);

Note that the initial condition should also be periodic on the domain.

using Plots
u_init = (x, y) -> sin(4 * π * x) * exp(cos(π * y))
U₀ = mtx(u_init)
M = maximum(abs, U₀)
contour(x, y, U₀';
    color=:redsblues,  clims=(-M, M), 
    aspect_ratio=1,
    xaxis=("x", (-1, 1)),  yaxis=("y", (-1, 1)), 
    title="Initial condition" )

This function computes the time derivative for the unknowns. The actual calculations take place using the matrix shape.

function du_dt(u, α, t)
    U = unvec(u)
    Uxx = Dxx * U
    Uyy = U * Dyy'            # 2nd partials
    du_dt = α * (Uxx + Uyy)    # PDE
    return vec(du_dt)
end;

Since this problem is parabolic, a stiff integrator is appropriate.

using OrdinaryDiffEq
IVP = ODEProblem(du_dt, vec(U₀), (0, 0.2), 0.1)
sol = solve(IVP, Rodas4P());

Here is an animation of the solution.

anim = @animate for t in range(0, 0.2, 81)
    surface(x, y, unvec(sol(t))';
        color=:redsblues,  clims=(-M, M),
        xaxis=(L"x", (-1, 1)), 
        yaxis=(L"y", (-1, 1)), 
        zlims=(-M, M),
        title=@sprintf("Heat equation, t=%.3f", t),
        dpi=150, colorbar=:none)
end
mp4(anim, "figures/2d-heat.mp4");

Example 13.2.3

The first step is to define a discretization of the domain.

m, n = 50, 36
x, Dx, Dxx = FNC.diffcheb(m, [-1, 1])
y, Dy, Dyy = FNC.diffcheb(n, [-1, 1])
mtx, X, Y, _ = FNC.tensorgrid(x, y)
U₀ = mtx( (x, y) -> (1 + y) * (1 - x)^4 * (1 + x)^2 * (1 - y^4) );

There are really two grids now: the full grid and the subset grid of interior points. Since the IVP unknowns are on the interior grid, that is the one we need to change shapes on. We also need the functions extend and chop to add and remove boundary values.

chop = U -> U[2:m, 2:n]
extend = U -> [zeros(m+1) [zeros(1, n-1); U; zeros(1, n-1)] zeros(m+1)]
unvec = u -> reshape(u, m-1, n-1)
pack = U -> vec(chop(U))
unpack = w -> extend(unvec(w))

#43 (generic function with 1 method)

Now we can define and solve the IVP using a stiff solver.

function dw_dt(w, ϵ, t)
    U = unpack(w)
    Ux, Uxx = Dx * U, Dxx * U
    Uyy = U * Dyy'
    du_dt = @. 1 - Ux + ϵ * (Uxx + Uyy)
    return pack(du_dt)
end

IVP = ODEProblem(dw_dt, pack(U₀), (0.0, 2), 0.05)
w = solve(IVP, Rodas4P());

When we evaluate the solution at a particular value of $t$ , we get a vector of the interior grid values. The same unpack function above converts this to a complete matrix of grid values.

U = t -> unpack(w(t))
contour(x, y, U(0.5)';
    fill=true,  color=:blues,  levels=20, l=0,
    aspect_ratio=1,  xlabel=L"x",  ylabel=L"y",
    title="Solution at t = 0.5")

anim = @animate for t in 0:0.02:2
    U = unpack(w(t))
    surface(x, y, U';
        layout=(1, 2),  size=(640, 320),
        xlabel=L"x",  ylabel=L"y",  zaxis=((0, 2), L"u(x,y)"),
        color=:blues,  alpha=0.66,  clims=(0, 2), colorbar=:none,
        title="Advection-diffusion",  dpi=150)
    contour!(x, y, U'; 
        levels=24, 
        aspect_ratio=1,  subplot=2, 
        xlabel=L"x",  ylabel=L"y",
        color=:blues,  clims=(0, 2),  colorbar=:none,
        title=@sprintf("t = %.2f", t))
end
mp4(anim, "figures/2d-advdiff.mp4");

Example 13.2.4

We start with the discretization and initial condition.

m, n = 40, 40
x, Dx, Dxx = FNC.diffcheb(m, [-2, 2])
y, Dy, Dyy = FNC.diffcheb(n, [-2, 2])
mtx, X, Y, _ = FNC.tensorgrid(x, y)
U₀ = mtx( (x, y) -> (x + 0.2) * exp(-12 * (x^2 + y^2)) )
V₀ = zeros(size(U₀));

Note that because $u$ is known on the boundary, while $v$ is unknown over the full grid, there are two different sizes of vec/unvec operations. We also need to define functions to pack grid unknowns into a vector and to unpack them. When the unknowns for $u$ are packed, the boundary values are chopped off, and these are restored when unpacking.

_, _, _, unvec_v, _ = FNC.tensorgrid(x, y)
_, _, _, unvec_u, _ = FNC.tensorgrid(x[2:m], y[2:n])
chop = U -> U[2:m, 2:n]
extend = U -> [zeros(m+1) [zeros(1, n-1); U; zeros(1, n-1)] zeros(m+1)]
pack = (U, V) -> [vec(chop(U)); vec(V)]
N = (m-1) * (n-1)    # number of interior unknowns
unpack = w -> ( extend(unvec_u(w[1:N])), unvec_v(w[N+1:end]) )

#55 (generic function with 1 method)

We can now define and solve the IVP. Since this problem is hyperbolic, not parabolic, a nonstiff integrator is faster than a stiff one.

function dw_dt(w, c, t)
    U, V = unpack(w)
    du_dt = V
    dv_dt = c^2 * (Dxx * U + U * Dyy')
    return pack(du_dt, dv_dt)
end

IVP = ODEProblem(dw_dt, pack(U₀, V₀), (0, 4.0), 1)
sol = solve(IVP, Tsit5())
U = t -> unpack(sol(t))[1]

#57 (generic function with 1 method)

anim = @animate for t in 0:4/100:4
    Ut = U(t)
    surface(x, y, Ut';
        layout=(1, 2), size=(640, 320),
        xlabel=L"x",  ylabel=L"y",  zaxis=((-0.1, 0.1), L"u(x,y)"),
        color=:redsblues,  alpha=0.66,  clims=(-0.1, 0.1), colorbar=:none,
        title="Wave equation",  dpi=150)
    contour!(x, y, Ut'; 
        levels=24,  subplot=2, 
        aspect_ratio=1,
        xlabel=L"x",  ylabel=L"y",
        color=:redsblues,  clims=(-0.1, 0.1), 
        colorbar=:none,  title=@sprintf("t = %.2f", t))
end
mp4(anim, "figures/2d-wave.mp4");

13.3 Laplace and Poisson equations¶

Example 13.3.1

A = [1 2; -2 0]

2×2 Matrix{Int64}:
  1  2
 -2  0

B = [1 10 100; -5 5 3]

2×3 Matrix{Int64}:
  1  10  100
 -5   5    3

Applying the definition manually, we get

A_kron_B = [
    A[1, 1]*B A[1, 2]*B;
    A[2, 1]*B A[2, 2]*B
    ]

4×6 Matrix{Int64}:
  1   10   100    2  20  200
 -5    5     3  -10  10    6
 -2  -20  -200    0   0    0
 10  -10    -6    0   0    0

That result should be the same as the following.

kron(A, B)

4×6 Matrix{Int64}:
  1   10   100    2  20  200
 -5    5     3  -10  10    6
 -2  -20  -200    0   0    0
 10  -10    -6    0   0    0

Example 13.3.2

We make a crude discretization for illustrative purposes.

m, n = 6, 5
x, Dx, Dxx = FNC.diffmat2(m, [0, 3])
y, Dy, Dyy = FNC.diffmat2(n, [-1, 1])
mtx, X, Y, unvec, is_boundary = FNC.tensorgrid(x, y)

(FNCFunctions.var"#mtx#93"{Vector{Float64}, Vector{Float64}}([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], [-1.0, -0.6, -0.19999999999999996, 0.20000000000000018, 0.6000000000000001, 1.0]), [0.0 0.0 … 0.0 0.0; 0.5 0.5 … 0.5 0.5; … ; 2.5 2.5 … 2.5 2.5; 3.0 3.0 … 3.0 3.0], [-1.0 -0.6 … 0.6000000000000001 1.0; -1.0 -0.6 … 0.6000000000000001 1.0; … ; -1.0 -0.6 … 0.6000000000000001 1.0; -1.0 -0.6 … 0.6000000000000001 1.0], FNCFunctions.var"#unvec#92"{Int64, Int64}(5, 6), Bool[1 1 … 1 1; 1 0 … 0 1; … ; 1 0 … 0 1; 1 1 … 1 1])

Next, we evaluate ϕ on the grid to get the forcing vector of the linear system.

ϕ = (x, y) -> x^2 - y + 2
b = vec(mtx(ϕ));

Here are the coefficients for the PDE collocation, before any modifications are made for the boundary conditions. The combination of Kronecker products and finite differences produces a characteristic sparsity pattern.

using SparseArrays, Plots
A = kron(I(n+1), sparse(Dxx)) + kron(sparse(Dyy), I(m+1))
spy(A;
    color=:blues,  m=3,
    title="System matrix before boundary conditions")

The number of equations is equal to $(m+1)(n+1)$ , which is the total number of points on the grid.

N = length(b)

42

The combination of Kronecker products and finite differences produces a characteristic sparsity pattern.

We now use the Boolean array that indicates where the boundary points lie in the grid.

spy(sparse(is_boundary);
    m=3,  color=:darkblue, 
    legend=:none,  title="Boundary points",
    xaxis=("column index", [0, n+2]), 
    yaxis=("row index", [0, m+2]))

In order to impose Dirichlet boundary conditions, we replace the boundary rows of the system by rows of the identity.

I_N = I(N)
idx = vec(is_boundary)
A[idx, :] .= I_N[idx, :];     # Dirichlet conditions

spy(A;
    color=:blues,  m=3,
    title="System matrix with boundary conditions")

Finally, we must replace the rows in the vector $\mathbf{b}$ by the boundary values being assigned to the boundary points. Here, we let the boundary values be zero everywhere.

b[idx] .= 0;                 # Dirichlet values

Now we can solve for $\mathbf{u}$ and reinterpret it as the matrix-shaped $\mathbf{U}$ , the solution on our grid.

u = A \ b
U = unvec(u)

7×6 Matrix{Float64}:
 0.0   0.0        0.0        0.0        0.0       0.0
 0.0  -0.549304  -0.758277  -0.712098  -0.453391  0.0
 0.0  -0.917873  -1.30273   -1.24377   -0.798473  0.0
 0.0  -1.21928   -1.74064   -1.67911   -1.09539   0.0
 0.0  -1.39869   -1.97682   -1.91786   -1.27929   0.0
 0.0  -1.21024   -1.65843   -1.61225   -1.11433   0.0
 0.0   0.0        0.0        0.0        0.0       0.0

Example 13.3.3

First we define the problem on $[0,1]\times[0,2]$ .

f = (x, y) -> -sin(3x * y - 4y) * (9y^2 + (3x - 4)^2);
g = (x, y) -> sin(3x * y - 4y);
xspan = [0, 1];
yspan = [0, 2];

Here is the finite-difference solution.

x, y, U = FNC.poissonfd(f, g, 40, xspan, 60, yspan);

surface(x, y, U';
    color=:viridis,
    title="Solution of Poisson's equation",
    xaxis=(L"x"),  yaxis=(L"y"),  zaxis=(L"u(x,y)"),
    right_margin=3Plots.mm,  camera=(70, 50))

The error is a smooth function of $x$ and $y$ . It must be zero on the boundary; otherwise, we have implemented boundary conditions incorrectly.

error = [g(x, y) for x in x, y in y] - U;
M = maximum(abs, error)
contour(x, y, error';
    levels=17, 
    clims=(-M, M), color=:redsblues, 
    colorbar=:bottom,  aspect_ratio=1,
    title="Error", 
    xaxis=(L"x"),  yaxis=(L"y"),
    right_margin=7Plots.mm)
plot!([0, 1, 1, 0, 0], [0, 0, 2, 2, 0], l=(2, :black))

13.4 Nonlinear elliptic PDEs¶

Example 13.4.2

All we need to define are ϕ from (13.4.2) for the PDE, and a trivial zero function for the boundary condition.

λ = 1.5
ϕ = (X, Y, U, Ux, Uxx, Uy, Uyy) -> @. Uxx + Uyy - λ / (U + 1)^2;
g = (x, y) -> 0;

Here is the solution for $m=15$ , $n=8$ .

u = FNC.elliptic(ϕ, g, 15, [0, 2.5], 8, [0, 1]);

using Plots
x = range(0, 2.5, 100)
y = range(0, 1, 50)
U = [u(x, y) for x in x, y in y]
contourf(x, y, U';
    color=:blues,  l=0,
    aspect_ratio=1,
    xlabel=L"x",  ylabel=L"y",  zlabel=L"u(x,y)",
    title="Deflection of a MEMS membrane",
    right_margin=3Plots.mm)

In the absence of an exact solution, how can we be confident that the solution is accurate? First, the Levenberg iteration converged without issuing a warning, so we should feel confident that the discrete equations were solved. We can check the boundary values easily. For example,

x_test = range(0, 2.5, 100)
norm([u(x, 0) - g(x, 0) for x in x_test], Inf)

2.4918568873705487e-23

Assuming that we encoded the PDE correctly, the remaining source error is truncation from the discretization. We can estimate that by refining the grid a bit and seeing how much the numerical solution changes.

x_test = range(0, 2.5, 6)
y_test = range(0, 1, 6)
mtx_test, _ = FNC.tensorgrid(x_test, y_test)
mtx_test(u)

6×6 Matrix{Float64}:
  0.0           2.38731e-25  -4.92724e-24  …   7.74378e-25   0.0
  2.91976e-24  -0.147964     -0.226459        -0.147964     -4.94098e-25
 -9.59403e-24  -0.195861     -0.305949        -0.195861      1.55867e-24
 -1.77355e-23  -0.195861     -0.305949        -0.195861      1.30691e-23
 -8.41765e-25  -0.147964     -0.226459        -0.147964     -4.47272e-24
  0.0           2.5961e-24   -1.38235e-24  …   4.40206e-24   0.0

u = FNC.elliptic(ϕ, g, 25, [0, 2.5], 14, [0, 1]);
mtx_test(u)

6×6 Matrix{Float64}:
  0.0          -2.32202e-22   3.64848e-23  …  -7.02464e-24   0.0
 -9.25622e-23  -0.147958     -0.226453        -0.147958     -1.59215e-22
 -1.02101e-22  -0.195861     -0.305929        -0.195861     -1.96135e-22
 -1.47369e-21  -0.195861     -0.305929        -0.195861     -4.03143e-22
 -1.46345e-22  -0.147958     -0.226453        -0.147958      1.00683e-22
  0.0          -4.08671e-23   3.06048e-22  …  -4.38852e-23   0.0

The original solution seems to be accurate to about four digits.

Example 13.4.3

ϕ = (X, Y, U, Ux, Uxx, Uy, Uyy) -> @. 1 - Ux - 2Uy + 0.05 * (Uxx + Uyy)
g = (x, y) -> 0
u = FNC.elliptic(ϕ, g, 32, [-1, 1], 32, [-1, 1]);

x = y = range(-1, 1, 80)
U = [u(x, y) for x in x, y in y]
contourf(x, y, U';
    color=:viridis, 
    aspect_ratio=1,
    xlabel=L"x",  ylabel=L"y",  zlabel=L"u(x,y)",
    title="Steady advection–diffusion")

Example 13.4.4

The following defines the PDE and a nontrivial Dirichlet boundary condition for the square $[0,1]^2$ .

ϕ = (X, Y, U, Ux, Uxx, Uy, Uyy) -> @. U * (1 - U^2) + 0.05 * (Uxx + Uyy)
g = (x, y) -> tanh(5 * (x + 2y - 1));

We solve the PDE and then plot the result.

u = FNC.elliptic(ϕ, g, 36, [0, 1], 36, [0, 1]);

x = y = range(0, 1, 80)
U = [u(x, y) for x in x, y in y]
contourf(x, y, U';
    color=:viridis, 
    aspect_ratio=1,
    xlabel=L"x",  ylabel=L"y",  zlabel=L"u(x,y)", 
    title="Steady Allen-Cahn",
    right_margin=3Plots.mm)