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.

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

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, n = 80, 70;
x, Dx, _ = FNC.diffmat2(m, [0, 2]);
y, Dy, _ = FNC.diffmat2(n, [1, 3]);
mtx = (f, x, y) -> [f(x, y) for x in x, y in y];

Here is how the exact first partial derivatives look on the grid.

plot(layout=(1, 2), aspect_ratio=1, xlabel="x", ylabel="y")
contour!(x, y, mtx(∂u_∂x, x, y)';
    levels=12, colormap=:viridis, subplot=1, title="∂u/∂x")
contour!(x, y, mtx(∂u_∂y, x, y)';
    levels=12, colormap=:viridis, subplot=2, title="∂u/∂y")

Next, we compute the finite-difference approximations of these derivatives on the grid.

U = mtx(u, x, y)
∂xU = Dx * U
∂yU = U * Dy';

Finally, we compare the finite-difference approximations to the exact values, using a colormap that is symmetric about zero.

err_x = mtx(∂u_∂x, x, y) - ∂xU
err_y = mtx(∂u_∂y, x, y) - ∂yU

plot(layout=(1, 2), aspect_ratio=1, 
    xlabel="x", ylabel="y")
M = maximum(abs, err_x)
contour!(x, y, err_x';
    levels=15, clims=(-M, M), colormap=:redsblues, fill=true,
    l=false, subplot=1, title="∂u/∂x error")
M = maximum(abs, err_y)
contour!(x, y, err_y';
    levels=15, clims=(-M, M), colormap=:redsblues, fill=true,
    l=false, subplot=2, title="∂u/∂y error")

Not surprisingly, the errors are largest where the derivatives themselves are largest.

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}:
 7
 3
 1
 6
 7
 8

The unvec operation is the inverse of vec.

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

2×3 Matrix{Int64}:
 7  1  7
 3  6  8

Example 13.2.2

We start by defining the discretization of the rectangle.

m, n = (60, 25)
x, Dx, Dxx = FNC.diffper(m, [-1, 1])
y, Dy, Dyy = FNC.diffper(n, [-1, 1])
mtx, X, Y, unvec, _ = FNC.tensorgrid(x, y);

Here is the initial condition, evaluated on the grid.

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

The following function computes the time derivative for the unknowns, which have a vector shape. The actual calculations, however, 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 time integrator is appropriate.

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

We can use the function sol defined above to get the solution at any time. Its output is the matrix $\mathbf{U}$ of values on the grid. An animation shows convergence of the solution toward a uniform value.

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

We define functions extend and chop to deal with the Dirichlet boundary conditions.

chop = U -> U[2:m, 2:n]
z = zeros(1, n-1)
extend = W -> [zeros(m+1) [z; W; z] zeros(m+1)]

#31 (generic function with 1 method)

Next, we define the pack and unpack functions, using another call to Algorithm 13.2.1 to get reshaping functions for the interior points.

_, _, _, unvec, _ = FNC.tensorgrid(x[2:m], y[2:n])
pack = U -> vec(chop(U))
unpack = w -> extend(unvec(w))

#35 (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.3), L"u(x,y)"),
        color=:blues,  alpha=0.66,  clims=(0, 2.3), 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.3),  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₀));

We need to define chopping and extension for the $u$ component. This looks the same as in Example 13.2.3.

chop = U -> U[2:m, 2:n]
extend = U -> [zeros(m+1) [zeros(1, n-1); U; zeros(1, n-1)] zeros(m+1)]

#43 (generic function with 1 method)

While vec is the same for both the interior and full grids, the unvec operation is defined differently for them.

_, _, _, unvec_v, _ = FNC.tensorgrid(x, y)
_, _, _, unvec_u, _ = FNC.tensorgrid(x[2:m], y[2:n])
N = (m-1) * (n-1)    # number of interior unknowns

pack = (U, V) -> [vec(chop(U)); vec(V)]
unpack = w -> (extend(unvec_u(w[1:N])), unvec_v(w[N+1:end]));

We can now define and solve the IVP. Since this problem is hyperbolic, 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];

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

Here is a look at the matrix we called $\mathbf{L}$ (the discrete Laplacian), 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 = (m+1) * (n+1)

42

We now use the final output from Algorithm 13.2.1, which is a Boolean array indicating 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 use the boundary indicator to index into the rows of the system.

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

Next, we evaluate ϕ on the grid to get the forcing vector of the linear system, and then modify the boundary rows to hold the boundary values—in this case, zero.

ϕ = (x, y) -> x^2 - y + 2
b = vec(mtx(ϕ));
b[idx] .= 0;    # Dirichlet values

Now we can solve the linear system $\mathbf{A}\mathbf{u}=\mathbf{b}$ 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)