Nonlinear elliptic PDEs - Fundamentals of Numerical Computation

Many nonlinear elliptic PDEs include references to the Laplacian operator.

More generally, we want to solve the nonlinear PDE

\phi(x,y,u,u_x,u_y,u_{xx},u_{yy}) = 0

(13.4.2)

in the interior of a rectangle $R$ , subject to the Dirichlet condition

u(x,y) = g(x,y)

(13.4.3)

on the boundary of $R$ .

13.4.1Implementation¶

In order to solve for as few unknowns as possible, we use a Chebyshev discretization of the domain. The core idea is to formulate collocation equations at the grid points based on discrete approximations of (13.4.2) and (13.4.3). If the PDE is nonlinear, then these equations are also nonlinear and are to be solved by a quasi-Newton iteration. Algorithm 13.4.1 is our implementation.

Algorithm 13.4.1 (elliptic)

Julia

MATLAB

Python

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

Solution of elliptic PDE by Chebyshev collocation

elliptic.m

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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
function u = elliptic(phi, g, m, xspan, n, yspan)
%ELLIPTIC   Solve an elliptic PDE in 2d.
% Input:
%   phi          defines phi(x, y, u, u_x, u_xx, u_y, u_yy) = 0 (function)
%   g            boundary condition (function)
%   m, xspan     size and interval of x discretization (integer, 2-vector)
%   n, yspan     size and interval of y discretization (integer, 2-vector)
% Output:
%   U            solution (n+1 by n+1)
%   X,Y          coordinate matrices (n+1 by n+1)

    % Discretize the domain.
    [x, Dx, Dxx] = diffcheb(m, xspan);
    [y, Dy, Dyy] = diffcheb(n, yspan);
    [mtx, X, Y, vec, unvec, is_boundary] = tensorgrid(x, y);

    % Identify boundary locations and evaluate the boundary condition.
    idx = vec(is_boundary);
    gb = g(X(idx), Y(idx));

    % Evaluate the PDE+BC residual.
    function r = residual(u)
        U = unvec(u);
        R = phi(X, Y, U, Dx * U, Dxx * U, U * Dy', U * Dyy');  % PDE
        R(idx) = u(idx) - gb;                                  % boundary
        r = vec(R);
    end

    % Solve the equation.
    u = levenberg(@residual, vec(zeros(size(X))));
    U = unvec(u(:, end));

    function u = evaluate(xi, eta)
        v = zeros(1, n+1);
        for j = 1:n+1
            v(j) = chebinterp(x, U(:, j), xi);
        end
        u = chebinterp(y, v, eta);
    end

    u = @evaluate;
end

function f = chebinterp(x, v, xi)
    n = length(x) - 1;
    w = (-1.0) .^ (0:n)';
    w([1, n+1]) = w([1, n+1]) / 2;

    terms = w ./ (xi - x(:));
    if any(isinf(terms))     % exactly at a node
        f = v(xi == x);
    else
        f = sum(v(:) .* terms) / sum(terms);
    end
end

Solution of elliptic PDE by Chebyshev collocation

elliptic.py

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
def elliptic(f, g, m, xspan, n, yspan):
    """
    newtonpde(f, g, m, xspan, n, yspan)

    Newton's method with finite differences to solve the PDE f(u,x,y,disc)=0 on the
    rectangle xspan \times yspan, subject to g(x,y)=0 on the boundary. Use m+1
    points in x by n+1 points in y.

    Return matrices of the solution values, and the coordinate functions, on the grid.
    """
    from scipy.sparse.linalg import spsolve
    x, Dx, Dxx = diffcheb(m, xspan)
    y, Dy, Dyy = diffcheb(n, yspan)
    mtx, X, Y, vec, unvec, is_boundary = tensorgrid(x, y)

    # Evaluate the boundary condition at the boundary nodes.
    idx = vec(is_boundary)
    X_bd, Y_bd = vec(X)[idx], vec(Y)[idx]
    g_bd = g(X_bd, Y_bd)

    # Evaluate the PDE+BC residual.
    def residual(u):
        U = unvec(u)
        R = f(X, Y, U, Dx @ U, Dxx @ U, U @ Dy.T, U @ Dyy.T)    # PDE
        r = vec(R)
        r[idx] = u[idx] - g_bd                                  # BC
        return r
    
    # Solve the equation.
    u = levenberg(residual, vec(np.zeros(X.shape)))[-1]
    U = unvec(u)

    def evaluate(xi, eta):
        v = [chebinterp(y, u, eta) for u in U]
        return chebinterp(x, v, xi)
    
    return np.vectorize(evaluate)

Algorithm 13.4.1 first defines the discretization and then computes all the values of $g$ at the boundary nodes. It uses Algorithm 4.6.3 as the nonlinear solver, and it translates back and forth between vector and grid shapes for the unknowns. After the discrete PDE is collocated at the grid points, the boundary terms are replaced by the boundary residual.

Lines 38–41, which produce the value returned by Algorithm 13.4.1, provide a function that evaluates the numerical solution anywhere in the domain, as is explained next.

13.4.2Off-grid evaluation¶

A Chebyshev grid is clustered close to the boundary of the domain. As a result, simple piecewise interpolation to evaluate off the grid, as is done by plotting routines, may be unacceptably nonsmooth and inaccurate. Instead, we should use the global polynomial interpolation that is foundational to the Chebyshev spectral method.

Let $\mathbf{U}$ be a matrix of solution values on the Chebyshev grid, defining a function $u(x,y)$ , and let $(\xi,\eta)$ be a point where we want to evaluate $u(x,y)$ . Column $\mathbf{u}_j$ of the grid matrix represents values spanning all the $x_i$ while $y$ is fixed at $y_j$ . Therefore, we can define an interpolating polynomial $p_j(x)$ based on the values in $\mathbf{u}_j$ .

Now let $v_j = p_j(\xi)$ for $j=1,\ldots,n$ . The vector $\mathbf{v}$ is a discretization of $u(\xi,y)$ at the Chebyshev nodes in $y$ . It defines an interpolating polynomial $q(y)$ , and finally we have $u(\xi,\eta)=q(\eta)$ . You can think of the total process as reducing one dimension at a time through the action of evaluating a polynomial interpolant at a point.

The function returned by Algorithm 13.4.1 performs interpolation as described above, using a helper function chebinterp (not shown). The helper performs the evaluation of a polynomial interpolant in one variable using a modified implementation of Function 9.2.1 that exploits the barycentric weights for Chebyshev nodes given in (9.3.3).^[1]

Example 13.4.2 (MEMS model in 2D)

We solve the PDE

\Delta u - \frac{\lambda}{(u+1)^2} = 0

(13.4.4)

on the rectangle $[0,2.5] \times [0,1]$ , with a zero Dirichlet condition on the boundary.

Julia

MATLAB

Python

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)

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           1.88629e-22  -1.23684e-23  …  -1.24193e-22   0.0
 -1.32086e-25  -0.147958     -0.226453        -0.147958     -5.68645e-23
 -3.55135e-22  -0.195861     -0.305929        -0.195861      8.60294e-23
 -1.22576e-21  -0.195861     -0.305929        -0.195861     -2.26458e-22
 -2.49121e-22  -0.147958     -0.226453        -0.147958     -2.62429e-23
  0.0           1.52693e-22   2.54043e-22  …  -4.0251e-23    0.0

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

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.

lambda = 1.5;
phi = @(X, Y, U, Ux, Uxx, Uy, Uyy) Uxx + Uyy - lambda ./ (U + 1).^2;
g = @(x, y) zeros(size(x));

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

u = elliptic(phi, g, 15, [0, 2.5], 8, [0, 1]);
disp(sprintf("solution at (2, 0.6) is %.7f", u(2, 0.6)))

solution at (2, 0.6) is -0.2264594

Source

x = linspace(0, 2.5, 91);
y = linspace(0, 1, 51);
[mtx, X, Y] = tensorgrid(x, y);
clf,  pcolor(x, y, mtx(u)')
colormap(flipud(sky)),  shading interp,  colorbar
axis equal
xlabel("x"),  ylabel("y")
title("Deflection of a MEMS membrane")

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. Assuming that we encoded the PDE correctly, the remaining source of 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 = linspace(0, 2.5, 6);
y_test = linspace(0, 1 , 5);
mtx_test = tensorgrid(x_test, y_test);
format long
mtx_test(u)

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

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

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.

lamb = 1.5
phi = lambda x, y, u, ux, uxx, uy, uyy: uxx + uyy - lamb / (u + 1)**2
g = lambda x, y: 0

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

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

print(f"solution at (2, 0.6) is {u(2, 0.6):.7f}")

solution at (2, 0.6) is -0.2264594

Source

x = linspace(0, 2.5, 90)
y = linspace(0, 1, 60)
mtx, X, Y, _, _, _ = FNC.tensorgrid(x, y)
U = mtx(u)

pcolormesh(X.T, Y.T, U.T, cmap="viridis")
xlabel("$x$"),  ylabel("$y$"),  axis("equal")
colorbar()
title("Solution of the MEMS equation in 2d");

err = norm(u(x, 0) - g(x, 0), inf)
print(f"max error on bottom edge: {err:.2e}")

max error on bottom edge: 9.56e-23

x_test = linspace(0, 2.5, 6)
y_test = linspace(0, 1, 6)
mtx_test, X_test, Y_test, _, _, _ = FNC.tensorgrid(x_test, y_test)

with printoptions(precision=7, suppress=True):
    print(mtx_test(u))

[[ 0.         0.         0.         0.         0.         0.       ]
 [ 0.        -0.1479641 -0.2264594 -0.2264594 -0.1479641  0.       ]
 [ 0.        -0.1958612 -0.3059492 -0.3059492 -0.1958612  0.       ]
 [ 0.        -0.1958612 -0.3059492 -0.3059492 -0.1958612  0.       ]
 [ 0.        -0.1479641 -0.2264594 -0.2264594 -0.1479641  0.       ]
 [-0.         0.         0.         0.         0.         0.       ]]

u = FNC.elliptic(phi, g, 25, [0, 2.5], 14, [0, 1])
with printoptions(precision=7, suppress=True):
    print(mtx_test(u))

[[ 0.        -0.        -0.        -0.         0.         0.       ]
 [-0.        -0.1479584 -0.226453  -0.226453  -0.1479584  0.       ]
 [ 0.        -0.195861  -0.3059293 -0.3059293 -0.195861   0.       ]
 [-0.        -0.195861  -0.3059293 -0.3059293 -0.195861  -0.       ]
 [-0.        -0.1479584 -0.226453  -0.226453  -0.1479584 -0.       ]
 [-0.        -0.        -0.         0.        -0.         0.       ]]

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

Example 13.4.3 (Steady advection-diffusion in 2D)

The steady-state limit of an advection-diffusion equation is

1 - u_x - 2u_y + \epsilon \, \Delta u = 0.

(13.4.5)

Here we solve it with a homogeneous Dirichlet condition on the square $[-1,1]^2$ .

Julia

MATLAB

Python

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

phi = @(X, Y, U, Ux, Uxx, Uy, Uyy) 1 - Ux - 2*Uy + 0.05*(Uxx + Uyy);
g = @(x, y) zeros(size(x));
u = elliptic(phi, g, 32, [-1, 1], 32, [-1, 1]);

Source

x = linspace(-1, 1, 80);
y = x;
mtx = tensorgrid(x, y);
clf,  pcolor(x, y, mtx(u)')
colormap(parula),  shading interp,  colorbar
axis equal,  xlabel("x"),  ylabel("y")
title("Steady advection–diffusion")

Example 13.4.3

phi = lambda x, y, u, ux, uxx, uy, uyy: 1 - ux - 2*uy + 0.05 * (uxx + uyy)
g = lambda x, y: 0
u = FNC.elliptic(phi, g, 32, [-1, 1], 32, [-1, 1])

x = y = linspace(-1, 1, 70)
mtx, X, Y, _, _, _ = FNC.tensorgrid(x, y)
U = mtx(u)

pcolormesh(X.T, Y.T, U.T, cmap="viridis")
xlabel("$x$"),  ylabel("$y$"),  axis("equal")
colorbar()
title("Steady advection–diffusion");

Example 13.4.4 (Steady Allen–Cahn in 2D)

The stationary Allen–Cahn equation in two dimensions is

u(1-u^2)+\epsilon \, \Delta u = 0.

(13.4.6)

Julia

MATLAB

Python

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)

Example 13.4.4

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

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

We solve the PDE and then plot the result.

u = elliptic(phi, g, 36, [0, 1], 36, [0, 1]);

Source

x = linspace(0, 1, 80);
y = x;
mtx = tensorgrid(x, y);
clf,  pcolor(x, y, mtx(u)')
colormap(parula),  shading interp,  colorbar
axis equal,  xlabel("x"),  ylabel("y")
title("Steady Allen–Cahn")

Example 13.4.4

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

phi = lambda x, y, u, ux, uxx, uy, uyy: u * (1 - u**2) + 0.05 * (uxx + uyy)
g = lambda x, y: tanh(5 * (x + 2*y - 1))

We solve the PDE and then plot the result.

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

Source

x = y = linspace(0, 1, 70)
mtx, X, Y, _, _, _ = FNC.tensorgrid(x, y)
U = mtx(u)
pcolormesh(X.T, Y.T, U.T, cmap="viridis")
xlabel("$x$"),  ylabel("$y$"),  axis("equal")
colorbar(),  title("Steady Allen–Cahn equation");

13.4.3Exercises¶

Exercise 13.4.2

⌨ A soap film stretched on a wire frame above the $(x,y)$ plane assumes a shape $u(x,y)$ of minimum area and is governed by

\begin{align*} \operatorname{div} \, \left( \frac{\operatorname{grad} u}{\sqrt{1 + u_x^2 + u_y^2}} \right) &= 0 \text{ in region $R$},\\ u(x,y) &= g(x,y) \text{ on the boundary of $R$}. \end{align*}

(13.4.8)

Solve the equation on $[-1,1]^2$ with boundary value $u(x,y)=\tanh(y-2x)$ , and make a surface plot of the result. (Hints: Don’t try to rewrite the PDE. Instead, modify Algorithm 13.4.1 so that ϕ is called with arguments (U,Dx,Dy), and compute the PDE in the form given. Also, since convergence is difficult in this problem, use the boundary data over the whole domain as the initial value for Algorithm 4.6.3.)

Exercise 13.4.3

⌨ Modify Algorithm 13.4.1 to solve (13.4.2) on $[a,b] \times [c,d]$ with the mixed boundary conditions

u = 0, \text{ if } x=a \text{ or } y = d, \qquad \frac{\partial u}{\partial n} = 0, \text{ if } x=b \text{ or } y = c,

(13.4.9)

where $\frac{\partial}{\partial n}$ is the derivative in the direction of the outward normal. Either condition can be used at a corner point. (Hint: Define index vectors for each side of the boundary square.) Apply your solver to the PDE $\Delta u + \sin(3\pi x) = 0$ on $[0,1]^2$ , and make a contour plot of the solution. Why do the level curves intersect two of the sides only at right angles?

Footnotes¶

The interpolation algorithm in Function 13.4.1 is inefficient when $u$ is to be evaluated on a finer grid, as for plotting. A more careful version could re-use the same values $v_j = p_j(\xi)$ for multiple values of η.
↩

References¶

Pelesko, J. A., & Driscoll, T. A. (2006). The Effect of the Small-Aspect-Ratio Approximation on Canonical Electrostatic MEMS Models. Journal of Engineering Mathematics, 53(3–4), 239–252. 10.1007/s10665-005-9013-2

Preface

Laplace and Poisson equations

Preface

Next steps