Note: This page has been translated by MathWorks. Please click here

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

Before solving a nonlinear elliptic PDE, from the **Solve** menu
in the PDE app, select **Parameters**. Then,
select the **Use nonlinear solver** check box and
click **OK**. At the command line, use `solvepde`

.

The basic idea is to use Gauss-Newton iterations to solve the nonlinear equations. Say you are trying to solve the equation

* r*(

In the FEM setting you solve the weak form of * r*(

$$u(x)={\displaystyle \sum {U}_{j}}{\varphi}_{j}$$

where **x** represents a 2-D or
3-D point. Then multiply the equation by an arbitrary test function * ϕ_{i}*,
integrate on the domain Ω, and use Green's formula and the boundary
conditions to obtain

$$\begin{array}{l}0=\rho \left(U\right)={\displaystyle \sum _{j}({\displaystyle \underset{\Omega}{\int}\left(c\left(x,U\right)\nabla {\varphi}_{j}(x)\right)\text{\hspace{0.17em}}\cdot \text{\hspace{0.17em}}\nabla {\varphi}_{j}(x)+a\left(x,U\right){\varphi}_{j}(x){\varphi}_{i}(x)\text{\hspace{0.17em}}dx}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+{\displaystyle \underset{\partial \Omega}{\int}q\left(x,U\right){\varphi}_{j}(x){\varphi}_{i}(x)\text{\hspace{0.17em}}ds}){U}_{j}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\displaystyle \underset{\Omega}{\int}f\left(x,U\right){\varphi}_{i}(x)\text{\hspace{0.17em}}dx}-{\displaystyle \underset{\partial \Omega}{\int}g\left(x,U\right){\varphi}_{i}(x)\text{\hspace{0.17em}}ds}\end{array}$$

which has to hold for all indices * i*.

The residual vector * ρ*(

* ρ*(

where the matrices * K*,

–∇ · (* c*(

Assume that you have a guess *U*^{(n)} of
the solution. If *U*^{(n)} is
close enough to the exact solution, an improved approximation *U*^{(n+1)} is
obtained by solving the linearized problem

$$\frac{\partial \rho \left({U}^{(n)}\right)}{\partial U}\left({U}^{(n+1)}-{U}^{(n)}\right)=-\alpha \rho \left({U}^{(n)}\right)$$

where $$\alpha $$ is a positive number. (It is
not necessary that * ρ(U) = 0* have a solution
even if

It is well known that for sufficiently small $$\alpha $$

$$\Vert \rho \left({U}^{(n+1)}\right)\Vert <\Vert \rho \left({U}^{(n)}\right)\Vert $$

and

$${p}_{n}={\left(\frac{\partial \rho \left({U}^{(n)}\right)}{\partial U}\right)}^{-1}\rho \left({U}^{(n)}\right)$$

is called a descent direction for $$\Vert \rho (U)\Vert $$, where $$\Vert \cdot \Vert $$ is the *L*_{2}-norm.
The iteration is

$${U}^{(n+1)}={U}^{(n)}+\alpha {p}_{n},$$

where $$\alpha $$ ≤ 1 is chosen as large as possible such that the step has a reasonable descent.

The *Gauss-Newton method* is local, and convergence
is assured only when *U*^{(0)} is
close enough to the solution. In general, the first guess may be outside
the region of convergence. To improve convergence from bad initial
guesses, a *damping* strategy is implemented for
choosing α, the *Armijo-Goldstein line search*.
It chooses the largest damping coefficient α out of the sequence
1, 1/2, 1/4, . . . such that the following inequality holds:

$$\Vert \rho \left({U}^{(n)}\right)\Vert -\Vert \rho \left({U}^{(n)}\right)+\alpha {p}_{n}\Vert \ge \frac{\alpha}{2}\Vert \rho \left({U}^{(n)}\right)\Vert $$

which guarantees a reduction of the residual norm by at least 1 – $$\alpha $$/2. Each step of the line-search algorithm requires an evaluation of the residual $$\rho \left({U}^{(n)}+\alpha {p}_{n}\right)$$.

An important point of this strategy is that when *U*^{(n)} approaches
the solution, then $$\alpha $$→1 and thus
the convergence rate increases. If there is a solution to * ρ(U)
= 0*, the scheme ultimately recovers the quadratic
convergence rate of the standard Newton iteration.

Closely related to the preceding problem is the choice of the
initial guess *U*^{(0)}.
By default, the solver sets *U*^{(0)} and
then assembles the FEM matrices * K* and

*U*^{(1)} = *K*^{–1}*F*

The damped Gauss-Newton iteration is then started with *U*^{(1)},
which should be a better guess than *U*^{(0)}.
If the boundary conditions do not depend on the solution * u*,
then

There are situations where * U^{(0)} =
0* makes no sense or convergence is impossible.

In some situations you may already have a good approximation and the nonlinear solver can be started with it, avoiding the slow convergence regime. This idea is used in the adaptive mesh generator. It computes a solution $$\tilde{U}$$ on a mesh, evaluates the error, and may refine certain triangles. The interpolant of $$\tilde{U}$$ is a very good starting guess for the solution on the refined mesh.

In general the exact Jacobian

$${J}_{n}=\frac{\partial \rho \left({U}^{(n)}\right)}{\partial U}$$

is not available. Approximation of * J_{n}* by
finite differences in the following way is expensive but feasible.
The

$$\frac{\rho \left({U}^{(n)}+\epsilon {\varphi}_{i}\right)-\rho \left({U}^{(n)}\right)}{\epsilon}$$

which implies the assembling of the FEM matrices for the triangles
containing grid point * i*. A very simple approximation
to

*U*^{(n+1)} = *K*^{–1}* F *.

This is equivalent to approximating the Jacobian with the stiffness
matrix. Indeed, since * ρ*(

$${U}^{(n+1)}={U}^{(n)}-{J}_{n}^{-1}\rho \left({U}^{(n)}\right)={U}^{(n)}-{K}^{-1}\left(K{U}^{(n)}-F\right)={K}^{-1}F$$

In many cases the convergence rate is slow, but the cost of each iteration is cheap.

The Partial Differential Equation Toolbox™ nonlinear solver
also provides for a compromise between the two extremes. To compute
the derivative of the mapping * U*→

$$\begin{array}{c}\frac{\partial {\left(KU\right)}_{i}}{\partial {U}_{j}}=\underset{\epsilon \to 0}{\mathrm{lim}}\frac{1}{\epsilon}{\displaystyle \sum _{l}({\displaystyle \underset{\Omega}{\int}c\left(U+\epsilon {\varphi}_{j}\right)\nabla {\varphi}_{l}\nabla {\varphi}_{i}\text{\hspace{0.17em}}dx\left({U}_{l}+\epsilon {\delta}_{l,j}\right)}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}-{\displaystyle \underset{\Omega}{\int}c\left(U\right)\nabla {\varphi}_{l}\nabla {\varphi}_{i}\text{\hspace{0.17em}}dx{U}_{l}})\\ ={\displaystyle \underset{\Omega}{\int}c\left(U\right)\nabla {\varphi}_{j}\nabla {\varphi}_{i}\text{\hspace{0.17em}}dx}+{\displaystyle \sum _{l}{\displaystyle \underset{\Omega}{\int}{\varphi}_{j}\frac{\partial c}{\partial u}\nabla {\varphi}_{l}\nabla {\varphi}_{i}\text{\hspace{0.17em}}dx{U}_{l}}}\end{array}$$

The first integral term is nothing more than * K_{i,j}*.

The second term is "lumped," i.e., replaced by
a diagonal matrix that contains the row sums. Since Σ_{j}* ϕ_{j}* =
1, the second term is approximated by

$${\delta}_{i,j}{\displaystyle \sum _{l}{\displaystyle \underset{\Omega}{\int}\frac{\partial c}{\partial u}\nabla {\varphi}_{l}\nabla {\varphi}_{i}\text{\hspace{0.17em}}dx{U}_{l}}}$$

which is the * i*th component of

$$-{\displaystyle \underset{\Omega}{\int}\frac{\partial f}{\partial u}{\varphi}_{i}{\varphi}_{j}\text{\hspace{0.17em}}dx}$$

which is the mass matrix associated with the coefficient ∂* f*/∂

$$J={K}^{(c)}+{M}^{(a-{f}^{\prime})}+\text{diag}\left(\left({K}^{({c}^{\prime})}+{M}^{({a}^{\prime})}\right)U\right)$$

where the differentiation is with respect to * u*,

$$\begin{array}{l}-\nabla \cdot (c\nabla u)+(a-f\text{'})u=0\\ -\nabla \cdot (c\text{'}\nabla u)+a\text{'}u=0\end{array}$$

and then produces the approximate Jacobian. The differentiations of the coefficients are done numerically.

In the general setting of elliptic systems, the boundary conditions are appended to the stiffness matrix to form the full linear system:

$$\tilde{K}\tilde{U}=\left[\begin{array}{cc}K& {H}^{\prime}\\ H& 0\end{array}\right]\left[\begin{array}{c}U\\ \mu \end{array}\right]=\left[\begin{array}{c}F\\ R\end{array}\right]=\tilde{F}$$

where the coefficients of $$\tilde{K}$$ and $$\tilde{F}$$ may depend on the solution $$\tilde{U}$$. The "lumped" approach approximates the derivative mapping of the residual by

$$\left[\begin{array}{cc}J& {H}^{\prime}\\ H& 0\end{array}\right]$$

The nonlinearities of the boundary conditions and the dependencies
of the coefficients on the derivatives of $$\tilde{U}$$ are not properly linearized
by this scheme. When such nonlinearities are strong, the scheme reduces
to the fix-point iteration and may converge slowly or not at all.
When the boundary conditions are linear, they do not affect the convergence
properties of the iteration schemes. In the Neumann case they are
invisible (* H* is an empty matrix) and in the Dirichlet
case they merely state that the residual is zero on the corresponding
boundary points.

Was this topic helpful?