from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

13.4.Nonlinear elliptic PDEs¶

Many nonlinear elliptic PDEs include references to the Laplacian operator.

More generally, we want to solve the nonlinear PDE

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

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=0,\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]

13.4.3Exercises¶