Nonlinear elliptic PDEs

Many nonlinear elliptic PDEs include references to the Laplacian operator.

More generally, we want to solve the nonlinear equation

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

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. Function 13.4.1 is our implementation.

Function 13.4.1 first defines the discretization and then computes all the values of $g$ at the boundary nodes. It uses Function 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 Function 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, and the grid values may be accurate even for modest grid sizes. As a result, simple piecewise interpolation to evaluate off the grid, as is done by plotting routines, may be unacceptably 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 Function 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]

13.4.3Exercises¶

1. ⌨ (a) Solve for the steady state of

   $$u_t = - u_y - x - 2 + \epsilon ( u_{xx} + u_{yy} )$$

   (13.4.7)

   for $\epsilon=1$ in $[-1,1]\times[-1,1]$, subject to a homogeneous Dirichlet boundary condition. Use $m=n=30$ points and plot the solution.

   (b) Repeat part (a) for $\epsilon=0.1$, which weakens the influence of diffusion relative to advection.

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 Function 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 levenberg.)

3. Modify Function 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 domain.) 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?