# Linearization at the differential equation level
<div id="nonlin:pdelevel"></div>

The attention is now turned to nonlinear partial differential
equations (PDEs) and application of the techniques explained above for
ODEs.  The model problem is a nonlinear diffusion equation for
$u(\x,t)$:

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:model:pde"></div>

$$
\begin{equation}
\frac{\partial u}{\partial t} = \nabla\cdot (\dfc(u)\nabla u) + f(u),\quad
\x\in\Omega,\ t\in (0,T],
\label{nonlin:pdelevel:model:pde} \tag{1}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:model:Neumann"></div>

$$
\begin{equation}  
-\dfc(u)\frac{\partial u}{\partial n} = g,\quad \x\in\partial\Omega_N,\ 
t\in (0,T],
\label{nonlin:pdelevel:model:Neumann} \tag{2}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:model:Dirichlet"></div>

$$
\begin{equation}  
u = u_0,\quad \x\in\partial\Omega_D,\ t\in (0,T]\thinspace .
\label{nonlin:pdelevel:model:Dirichlet} \tag{3}
\end{equation}
$$

In the present section, our aim is to discretize this problem in time
and then present techniques for linearizing the time-discrete PDE
problem "at the PDE level" such that we transform the nonlinear
stationary PDE problem at each time level into a sequence of linear
PDE problems, which can be solved using any method for linear
PDEs. This strategy avoids the solution of systems of nonlinear
algebraic equations.  In the section [1D stationary nonlinear differential equations](#nonlin:alglevel:1D) we shall take
the opposite (and more common) approach: discretize the nonlinear
problem in time and space first, and then solve the resulting
nonlinear algebraic equations at each time level by the methods of
the section [nonlin:systems:alg](#nonlin:systems:alg).  Very often, the two approaches are
mathematically identical, so there is no preference from a
computational efficiency point of view.  The details of the ideas
sketched above will hopefully become clear through the forthcoming
examples.


## Explicit time integration
<div id="nonlin:pdelevel:explicit"></div>

The nonlinearities in the PDE are trivial to deal with if we choose an
explicit time integration method for ([1](#nonlin:pdelevel:model:pde)),
such as the Forward Euler method:

$$
[D_t^+ u = \nabla\cdot (\dfc(u)\nabla u) + f(u)]^n,
$$

or written out,

$$
\frac{u^{n+1} - u^n}{\Delta t} = \nabla\cdot (\dfc(u^n)\nabla u^n)
+ f(u^n),
$$

which is a linear equation in the unknown $u^{n+1}$ with solution

$$
u^{n+1} = u^n + \Delta t\nabla\cdot (\dfc(u^n)\nabla u^n) +
\Delta t f(u^n)\thinspace .
$$

The disadvantage with this discretization is
the strict stability criterion $\Delta t \leq h^2/(6\max\alpha)$
for the case $f=0$ and a standard 2nd-order finite difference discretization
in 3D space with mesh cell sizes $h=\Delta x=\Delta y=\Delta z$.

<!-- BC -->

## Backward Euler scheme and Picard iteration
<div id="nonlin:pdelevel:Picard"></div>

A Backward Euler scheme for ([1](#nonlin:pdelevel:model:pde))
reads

$$
[D_t^- u = \nabla\cdot (\dfc(u)\nabla u) + f(u)]^n\thinspace .
$$

Written out,

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:pde:BE"></div>

$$
\begin{equation}
\frac{u^{n} - u^{n-1}}{\Delta t} = \nabla\cdot (\dfc(u^n)\nabla u^n)
+ f(u^n)\thinspace .
\label{nonlin:pdelevel:pde:BE} \tag{4}
\end{equation}
$$

This is a nonlinear PDE for the unknown function $u^n(\x)$. Such a
PDE can be viewed as a time-independent PDE where
$u^{n-1}(\x)$ is a known function.

We introduce a Picard iteration with $k$ as iteration counter.
A typical linearization of the $\nabla\cdot(\dfc(u^n)\nabla u^n)$ term
in iteration $k+1$ is to use the previously computed $u^{n,k}$
approximation in the diffusion coefficient: $\dfc(u^{n,k})$.
The nonlinear source term is treated similarly: $f(u^{n,k})$.
The unknown function $u^{n,k+1}$ then fulfills the linear PDE

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:pde:BE:Picard:k"></div>

$$
\begin{equation}
\frac{u^{n,k+1} - u^{n-1}}{\Delta t} = \nabla\cdot (\dfc(u^{n,k})
\nabla u^{n,k+1})
+ f(u^{n,k})\thinspace .
\label{nonlin:pdelevel:pde:BE:Picard:k} \tag{5}
\end{equation}
$$

The initial guess for the Picard iteration at this time level can be
taken as the solution at the previous time level: $u^{n,0}=u^{n-1}$.

We can alternatively apply the implementation-friendly
notation where $u$ corresponds to
the unknown we want to solve for, i.e., $u^{n,k+1}$ above, and $u^{-}$
is the most recently computed value, $u^{n,k}$ above. Moreover,
$u^{(1)}$ denotes the unknown function at the previous time level, $u^{n-1}$
above. The PDE to be solved in a Picard iteration then looks like

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:pde:BE:Picard"></div>

$$
\begin{equation}
\frac{u - u^{(1)}}{\Delta t} = \nabla\cdot (\dfc(u^{-})
\nabla u)
+ f(u^{-})\thinspace .
\label{nonlin:pdelevel:pde:BE:Picard} \tag{6}
\end{equation}
$$

At the beginning of the iteration we start with the value from the
previous time level: $u^{-}=u^{(1)}$, and after each
iteration, $u^{-}$ is updated to $u$.

**Remark on notation.**

The previous derivations of the numerical scheme for time discretizations
of PDEs have, strictly
speaking, a somewhat sloppy notation, but it is much used and convenient
to read. A more precise notation must
distinguish clearly between the exact solution of the PDE problem,
here denoted $\uex(\x,t)$, and the exact solution of the spatial
problem, arising after time discretization at each time level,
where ([4](#nonlin:pdelevel:pde:BE)) is an example. The latter
is here represented as $u^n(\x)$ and is an approximation to
$\uex(\x,t_n)$. Then we have another approximation $u^{n,k}(\x)$
to $u^n(\x)$ when solving the nonlinear PDE problem for
$u^n$ by iteration methods, as in ([5](#nonlin:pdelevel:pde:BE:Picard:k)).

In our notation, $u$ is a synonym for $u^{n,k+1}$ and $u^{(1)}$ is
a synonym for $u^{n-1}$, inspired by what are natural variable names
in a code.
We will usually state the PDE problem in terms of $u$ and
quickly redefine the symbol $u$ to mean the numerical approximation,
while $\uex$ is not explicitly introduced unless we need to talk about
the exact solution and the approximate solution at the same time.



## Backward Euler scheme and Newton's method
<div id="nonlin:pdelevel:Newton"></div>


At time level $n$, we have to solve the stationary PDE
([4](#nonlin:pdelevel:pde:BE)). In the previous section, we
saw how this can be done with Picard iterations.
Another alternative is to apply the idea of Newton's method
in a clever way.
Normally, Newton's method is defined for systems of *algebraic equations*,
but the idea of the method can be applied at the PDE level too.

### Linearization via Taylor expansions

Let $u^{n,k}$ be an approximation to the unknown $u^n$. We seek a
better approximation on
the form

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:Newton:ansatz"></div>

$$
\begin{equation}
u^{n} = u^{n,k} + \delta u\thinspace .
\label{nonlin:pdelevel:Newton:ansatz} \tag{7}
\end{equation}
$$

The idea is to insert ([7](#nonlin:pdelevel:Newton:ansatz)) in
([4](#nonlin:pdelevel:pde:BE)), Taylor expand the nonlinearities
and keep only the terms that are
linear in $\delta u$ (which makes ([7](#nonlin:pdelevel:Newton:ansatz))
an approximation for $u^{n}$). Then we can solve a linear PDE for
the correction $\delta u$ and use ([7](#nonlin:pdelevel:Newton:ansatz))
to find a new approximation

$$
u^{n,k+1}=u^{n,k}+\delta u
$$

to $u^{n}$.
Repeating this procedure gives a sequence $u^{n,k+1}$, $k=0,1,\ldots$
that hopefully converges to the goal $u^n$.

Let us carry out all the mathematical details for the nonlinear diffusion
PDE discretized by the Backward Euler method.
Inserting ([7](#nonlin:pdelevel:Newton:ansatz)) in
([4](#nonlin:pdelevel:pde:BE)) gives

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:pde:BE:Newton1"></div>

$$
\begin{equation}
\frac{u^{n,k} +\delta u - u^{n-1}}{\Delta t} =
\nabla\cdot (\dfc(u^{n,k} + \delta u)\nabla (u^{n,k}+\delta u))
+ f(u^{n,k}+\delta u)\thinspace .
\label{nonlin:pdelevel:pde:BE:Newton1} \tag{8}
\end{equation}
$$

We can Taylor expand $\dfc(u^{n,k} + \delta u)$ and
$f(u^{n,k}+\delta u)$:

$$
\begin{align*}
\dfc(u^{n,k} + \delta u) & = \dfc(u^{n,k}) + \frac{d\dfc}{du}(u^{n,k})
\delta u + \Oof{\delta u^2}\approx \dfc(u^{n,k}) + \dfc^{\prime}(u^{n,k})\delta u,\\ 
f(u^{n,k}+\delta u) &=  f(u^{n,k}) + \frac{df}{du}(u^{n,k})\delta u
+ \Oof{\delta u^2}\approx f(u^{n,k}) + f^{\prime}(u^{n,k})\delta u\thinspace .
\end{align*}
$$

Inserting the linear approximations of $\dfc$ and $f$ in
([8](#nonlin:pdelevel:pde:BE:Newton1)) results in

$$
\frac{u^{n,k} +\delta u - u^{n-1}}{\Delta t} =
\nabla\cdot (\dfc(u^{n,k})\nabla u^{n,k}) + f(u^{n,k}) + \nonumber
$$

$$
\qquad \nabla\cdot (\dfc(u^{n,k})\nabla \delta u)
+ \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla u^{n,k}) + \nonumber
$$

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:pde:BE:Newton2"></div>

$$
\begin{equation}  
\qquad \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla \delta u)
+ f^{\prime}(u^{n,k})\delta u\thinspace .
\label{nonlin:pdelevel:pde:BE:Newton2} \tag{9}
\end{equation}
$$

The term $\dfc^{\prime}(u^{n,k})\delta u\nabla \delta u$ is of
order $\delta u^2$
and therefore omitted since we expect the correction $\delta u$
to be small ($\delta u \gg \delta u^2$).
Reorganizing the equation gives a PDE
for $\delta u$ that we can write in short form as

$$
\delta F(\delta u; u^{n,k}) = -F(u^{n,k}),
$$

where

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:pde:BE:Newton2:F"></div>

$$
\begin{equation}
F(u^{n,k}) = \frac{u^{n,k} - u^{n-1}}{\Delta t} -
\nabla\cdot (\dfc(u^{n,k})\nabla u^{n,k}) + f(u^{n,k}),
\label{nonlin:pdelevel:pde:BE:Newton2:F} \tag{10}
\end{equation}
$$

$$
\delta F(\delta u; u^{n,k}) =
- \frac{1}{\Delta t}\delta u +
\nabla\cdot (\dfc(u^{n,k})\nabla \delta u) + \nonumber
$$

<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>

$$
\begin{equation}  
\quad \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla u^{n,k})
+ f^{\prime}(u^{n,k})\delta u\thinspace .
\label{_auto1} \tag{11}
\end{equation}
$$

Note that $\delta F$ is a linear function of $\delta u$, and
$F$ contains only terms that are known, such that
the PDE for $\delta u$ is indeed linear.

**Observations.**

The notational form $\delta F = -F$ resembles the Newton system $J\delta u =-F$
for systems of algebraic equations, with $\delta F$ as $J\delta u$.
The unknown vector in a linear system of algebraic equations enters
the system as a linear operator in terms of a
matrix-vector product ($J\delta u$), while at
the PDE level we have a linear differential operator instead
($\delta F$).



### Similarity with Picard iteration

We can rewrite the PDE for $\delta u$ in a slightly different way too
if we define $u^{n,k} + \delta u$ as $u^{n,k+1}$.

$$
\frac{u^{n,k+1} - u^{n-1}}{\Delta t} =
\nabla\cdot (\dfc(u^{n,k})\nabla u^{n,k+1}) + f(u^{n,k})\nonumber
$$

<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>

$$
\begin{equation}  
\qquad  + \nabla\cdot (\dfc^{\prime}(u^{n,k})\delta u\nabla u^{n,k})
+ f^{\prime}(u^{n,k})\delta u\thinspace .
\label{_auto2} \tag{12}
\end{equation}
$$

Note that the first line is the same PDE as arises in the Picard
iteration, while the remaining terms arise from the differentiations
that are an inherent ingredient in Newton's method.

### Implementation

For coding we want to introduce $u$ for $u^n$, $u^{-}$ for $u^{n,k}$ and
$u^{(1)}$ for $u^{n-1}$. The formulas for $F$ and $\delta F$
are then more clearly written as

<!-- Equation labels as ordinary links -->
<div id="nonlin:pdelevel:pde:BE:Newton2:F2"></div>

$$
\begin{equation}
F(u^{-}) = \frac{u^{-} - u^{(1)}}{\Delta t} -
\nabla\cdot (\dfc(u^{-})\nabla u^{-}) + f(u^{-}),
\label{nonlin:pdelevel:pde:BE:Newton2:F2} \tag{13}
\end{equation}
$$

$$
\delta F(\delta u; u^{-}) =
- \frac{1}{\Delta t}\delta u +
\nabla\cdot (\dfc(u^{-})\nabla \delta u) + \nonumber
$$

<!-- Equation labels as ordinary links -->
<div id="_auto3"></div>

$$
\begin{equation}  
\quad \nabla\cdot (\dfc^{\prime}(u^{-})\delta u\nabla u^{-})
+ f^{\prime}(u^{-})\delta u\thinspace .
\label{_auto3} \tag{14}
\end{equation}
$$

The form that orders the PDE as the Picard iteration terms plus
the Newton method's derivative terms becomes

$$
\frac{u - u^{(1)}}{\Delta t} =
\nabla\cdot (\dfc(u^{-})\nabla u) + f(u^{-}) + \nonumber
$$

<!-- Equation labels as ordinary links -->
<div id="_auto4"></div>

$$
\begin{equation}  
\qquad  \gamma(\nabla\cdot (\dfc^{\prime}(u^{-})(u - u^{-})\nabla u^{-})
+ f^{\prime}(u^{-})(u - u^{-}))\thinspace .
\label{_auto4} \tag{15}
\end{equation}
$$

The Picard and full Newton versions correspond to
$\gamma=0$ and $\gamma=1$, respectively.

### Derivation with alternative notation

Some may prefer to derive the linearized PDE for $\delta u$ using
the more compact notation. We start with inserting $u^n=u^{-}+\delta u$
to get

$$
\frac{u^{-} +\delta u - u^{n-1}}{\Delta t} =
\nabla\cdot (\dfc(u^{-} + \delta u)\nabla (u^{-}+\delta u))
+ f(u^{-}+\delta u)\thinspace .
$$

Taylor expanding,

$$
\begin{align*}
\dfc(u^{-} + \delta u) & \approx \dfc(u^{-}) + \dfc^{\prime}(u^{-})\delta u,\\ 
f(u^{-}+\delta u) & \approx f(u^{-}) + f^{\prime}(u^{-})\delta u,
\end{align*}
$$

and inserting these expressions gives a less cluttered PDE for $\delta u$:

$$
\begin{align*}
\frac{u^{-} +\delta u - u^{n-1}}{\Delta t} &=
\nabla\cdot (\dfc(u^{-})\nabla u^{-}) + f(u^{-}) + \\ 
&\qquad \nabla\cdot (\dfc(u^{-})\nabla \delta u)
+ \nabla\cdot (\dfc^{\prime}(u^{-})\delta u\nabla u^{-}) + \\ 
&\qquad \nabla\cdot (\dfc^{\prime}(u^{-})\delta u\nabla \delta u)
+ f^{\prime}(u^{-})\delta u\thinspace .
\end{align*}
$$

## Crank-Nicolson discretization
<div id="nonlin:pdelevel:Picard:CN"></div>

A Crank-Nicolson discretization of
([1](#nonlin:pdelevel:model:pde)) applies a centered difference
at $t_{n+\frac{1}{2}}$:

$$
[D_t u = \nabla\cdot (\dfc(u)\nabla u) + f(u)]^{n+\frac{1}{2}}\thinspace .
$$

The standard technique is to apply an arithmetic average for
quantities defined between two mesh points, e.g.,

$$
u^{n+\frac{1}{2}}\approx \frac{1}{2}(u^n + u^{n+1})\thinspace .
$$

However, with nonlinear terms we have many choices of formulating
an arithmetic mean:

<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>

$$
\begin{equation}
[f(u)]^{n+\frac{1}{2}} \approx f(\frac{1}{2}(u^n + u^{n+1}))
= [f(\overline{u}^t)]^{n+\frac{1}{2}},
\label{_auto5} \tag{16}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto6"></div>

$$
\begin{equation}  
[f(u)]^{n+\frac{1}{2}} \approx \frac{1}{2}(f(u^n) + f(u^{n+1}))
=[\overline{f(u)}^t]^{n+\frac{1}{2}},
\label{_auto6} \tag{17}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto7"></div>

$$
\begin{equation}  
[\dfc(u)\nabla u]^{n+\frac{1}{2}} \approx
\dfc(\frac{1}{2}(u^n + u^{n+1}))\nabla (\frac{1}{2}(u^n + u^{n+1}))
= [\dfc(\overline{u}^t)\nabla \overline{u}^t]^{n+\frac{1}{2}},
\label{_auto7} \tag{18}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto8"></div>

$$
\begin{equation}  
[\dfc(u)\nabla u]^{n+\frac{1}{2}} \approx
\frac{1}{2}(\dfc(u^n) + \dfc(u^{n+1}))\nabla (\frac{1}{2}(u^n + u^{n+1}))
= [\overline{\dfc(u)}^t\nabla\overline{u}^t]^{n+\frac{1}{2}},
\label{_auto8} \tag{19}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto9"></div>

$$
\begin{equation}  
[\dfc(u)\nabla u]^{n+\frac{1}{2}} \approx
\frac{1}{2}(\dfc(u^n)\nabla u^n + \dfc(u^{n+1})\nabla u^{n+1})
= [\overline{\dfc(u)\nabla u}^t]^{n+\frac{1}{2}}\thinspace .
\label{_auto9} \tag{20}
\end{equation}
$$

A big question is whether there are significant differences in accuracy
between taking the products of arithmetic means or taking the arithmetic
mean of products. [nonlin:exer:products:arith:mean](#nonlin:exer:products:arith:mean) investigates
this question, and the answer is that the approximation is
$\Oof{\Delta t^2}$ in both cases.


# 1D stationary nonlinear differential equations
<div id="nonlin:alglevel:1D"></div>

The section [Linearization at the differential equation level](#nonlin:pdelevel) presented methods for linearizing
time-discrete PDEs directly prior to discretization in space.  We can
alternatively carry out the discretization in space of the
time-discrete nonlinear PDE problem and get a system of nonlinear
algebraic equations, which can be solved by Picard iteration or
Newton's method as presented in the section [nonlin:systems:alg](#nonlin:systems:alg).
This latter approach will now be described in detail.

We shall work with the 1D problem

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:pde"></div>

$$
\begin{equation}
-(\dfc(u)u^{\prime})^{\prime} + au = f(u),\quad x\in (0,L),
\quad \dfc(u(0))u^{\prime}(0) = C,\ u(L)=D
\thinspace .
\label{nonlin:alglevel:1D:pde} \tag{21}
\end{equation}
$$

The problem ([21](#nonlin:alglevel:1D:pde)) arises from the stationary
limit of a diffusion equation,

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:pde:tver"></div>

$$
\begin{equation}
\frac{\partial u}{\partial t} = \frac{\partial}{\partial x}\left(
\alpha(u)\frac{\partial u}{\partial x}\right) - au + f(u),
\label{nonlin:alglevel:1D:pde:tver} \tag{22}
\end{equation}
$$

as $t\rightarrow\infty$ and $\partial u/\partial t\rightarrow 0$.
Alternatively, the problem ([21](#nonlin:alglevel:1D:pde)) arises
at each time level from implicit time discretization of
([22](#nonlin:alglevel:1D:pde:tver)). For example, a Backward Euler
scheme for ([22](#nonlin:alglevel:1D:pde:tver)) leads to

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:pde:tver:BE"></div>

$$
\begin{equation}
\frac{u^{n}-u^{n-1}}{\Delta t} =
\frac{d}{dx}\left(
\alpha(u^n)\frac{du^n}{dx}\right) - au^n + f(u^n)\thinspace .
\label{nonlin:alglevel:1D:pde:tver:BE} \tag{23}
\end{equation}
$$

Introducing $u(x)$ for $u^n(x)$, $u^{(1)}$ for $u^{n-1}$, and defining $f(u)$
in ([21](#nonlin:alglevel:1D:pde)) to be $f(u)$ in
([23](#nonlin:alglevel:1D:pde:tver:BE)) plus $u^{n-1}/\Delta t$, gives
([21](#nonlin:alglevel:1D:pde)) with $a=1/\Delta t$.


## Finite difference discretization
<div id="nonlin:alglevel:1D:fd"></div>


The nonlinearity in the differential equation
([21](#nonlin:alglevel:1D:pde)) poses no more difficulty than a variable
coefficient, as in the term $(\dfc(x)u^{\prime})^{\prime}$.  We can
therefore use a standard finite difference approach when discretizing
the Laplace term with a variable coefficient:

$$
[-D_x\dfc D_x u +au = f]_i\thinspace .
$$

Writing this out for a uniform mesh with points $x_i=i\Delta x$,
$i=0,\ldots,N_x$, leads to

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:deq0"></div>

$$
\begin{equation}
-\frac{1}{\Delta x^2}
\left(\dfc_{i+\frac{1}{2}}(u_{i+1}-u_i) -
\dfc_{i-\frac{1}{2}}(u_{i}-u_{i-1})\right)
+ au_i = f(u_i)\thinspace .
\label{nonlin:alglevel:1D:fd:deq0} \tag{24}
\end{equation}
$$

This equation is valid at all the mesh points $i=0,1,\ldots,N_x-1$.
At $i=N_x$ we have the Dirichlet condition $u_i=0$.
The only difference from the case with $(\dfc(x)u^{\prime})^{\prime}$ and $f(x)$ is that
now $\dfc$ and $f$ are functions of $u$ and not only of $x$:
$(\dfc(u(x))u^{\prime})^{\prime}$ and $f(u(x))$.

The quantity $\dfc_{i+\frac{1}{2}}$, evaluated between two mesh points,
needs a comment. Since $\dfc$ depends on $u$ and $u$ is only known
at the mesh points, we need to express $\dfc_{i+\frac{1}{2}}$ in
terms of $u_i$ and $u_{i+1}$. For this purpose we use an arithmetic
mean, although a harmonic mean is also common in this context if
$\dfc$ features large jumps.
There are two choices of arithmetic means:

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:dfc:mean:u"></div>

$$
\begin{equation}
\dfc_{i+\frac{1}{2}} \approx
\dfc(\frac{1}{2}(u_i + u_{i+1}) =
[\dfc(\overline{u}^x)]^{i+\frac{1}{2}},
\label{nonlin:alglevel:1D:fd:dfc:mean:u} \tag{25}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:dfc:mean:dfc"></div>

$$
\begin{equation}  
\dfc_{i+\frac{1}{2}} \approx
\frac{1}{2}(\dfc(u_i) + \dfc(u_{i+1})) = [\overline{\dfc(u)}^x]^{i+\frac{1}{2}}
\label{nonlin:alglevel:1D:fd:dfc:mean:dfc} \tag{26}
\end{equation}
$$

Equation ([24](#nonlin:alglevel:1D:fd:deq0)) with
the latter approximation then looks like

$$
-\frac{1}{2\Delta x^2}
\left((\dfc(u_i)+\dfc(u_{i+1}))(u_{i+1}-u_i) -
(\dfc(u_{i-1})+\dfc(u_{i}))(u_{i}-u_{i-1})\right)\nonumber
$$

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:deq"></div>

$$
\begin{equation}  
\qquad\qquad + au_i = f(u_i),
\label{nonlin:alglevel:1D:fd:deq} \tag{27}
\end{equation}
$$

or written more compactly,

$$
[-D_x\overline{\dfc}^x D_x u +au = f]_i\thinspace .
$$

At mesh point $i=0$ we have the boundary condition $\dfc(u)u^{\prime}=C$,
which is discretized by

$$
[\dfc(u)D_{2x}u = C]_0,
$$

meaning

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:Neumann:x0"></div>

$$
\begin{equation}
\dfc(u_0)\frac{u_{1} - u_{-1}}{2\Delta x} = C\thinspace .
\label{nonlin:alglevel:1D:fd:Neumann:x0} \tag{28}
\end{equation}
$$

The fictitious value $u_{-1}$ can be eliminated with the aid
of ([27](#nonlin:alglevel:1D:fd:deq)) for $i=0$.
Formally, ([27](#nonlin:alglevel:1D:fd:deq)) should be solved with
respect to $u_{i-1}$ and that value (for $i=0$) should be inserted in
([28](#nonlin:alglevel:1D:fd:Neumann:x0)), but it is algebraically
much easier to do it the other way around. Alternatively, one can
use a ghost cell $[-\Delta x,0]$ and update the $u_{-1}$ value
in the ghost cell according to ([28](#nonlin:alglevel:1D:fd:Neumann:x0))
after every Picard or Newton iteration. Such an approach means that
we use a known $u_{-1}$ value in ([27](#nonlin:alglevel:1D:fd:deq))
from the previous iteration.

## Solution of algebraic equations

### The structure of the equation system

The nonlinear algebraic equations ([27](#nonlin:alglevel:1D:fd:deq)) are
of the form $A(u)u = b(u)$ with

$$
\begin{align*}
A_{i,i} &= \frac{1}{2\Delta x^2}(\dfc(u_{i-1}) + 2\dfc(u_{i})
\dfc(u_{i+1})) + a,\\ 
A_{i,i-1} &= -\frac{1}{2\Delta x^2}(\dfc(u_{i-1}) + \dfc(u_{i})),\\ 
A_{i,i+1} &= -\frac{1}{2\Delta x^2}(\dfc(u_{i}) + \dfc(u_{i+1})),\\ 
b_i &= f(u_i)\thinspace .
\end{align*}
$$

The matrix $A(u)$ is tridiagonal: $A_{i,j}=0$ for $j > i+1$ and $j < i-1$.

The above expressions are valid for internal mesh points $1\leq i\leq N_x-1$.
For $i=0$ we need to express $u_{i-1}=u_{-1}$ in terms of $u_1$ using
([28](#nonlin:alglevel:1D:fd:Neumann:x0)):

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:Neumann:x0:um1"></div>

$$
\begin{equation}
u_{-1} = u_1 -\frac{2\Delta x}{\dfc(u_0)}C\thinspace .
\label{nonlin:alglevel:1D:fd:Neumann:x0:um1} \tag{29}
\end{equation}
$$

This value must be inserted in $A_{0,0}$. The expression for $A_{i,i+1}$
applies for $i=0$, and $A_{i,i-1}$ does not enter the system when $i=0$.

Regarding the last equation, its form depends on whether we include
the Dirichlet condition $u(L)=D$, meaning $u_{N_x}=D$, in the
nonlinear algebraic equation system or not. Suppose we choose
$(u_0,u_1,\ldots,u_{N_x-1})$ as unknowns, later referred to as
*systems without Dirichlet conditions*. The last equation
corresponds to $i=N_x-1$. mathcal{I}_t involves the boundary value $u_{N_x}$,
which is substituted by $D$. If the unknown vector includes the
boundary value, $(u_0,u_1,\ldots,u_{N_x})$, later referred to as
*system including Dirichlet conditions*, the equation for $i=N_x-1$
just involves the unknown $u_{N_x}$, and the final equation becomes
$u_{N_x}=D$, corresponding to $A_{i,i}=1$ and $b_i=D$ for $i=N_x$.

### Picard iteration

The obvious Picard iteration scheme is to use previously computed
values of $u_i$ in $A(u)$ and $b(u)$, as described more in detail in
the section [nonlin:systems:alg](#nonlin:systems:alg). With the notation $u^{-}$ for the
most recently computed value of $u$, we have the system $F(u)\approx
\hat F(u) = A(u^{-})u - b(u^{-})$, with $F=(F_0,F_1,\ldots,F_m)$,
$u=(u_0,u_1,\ldots,u_m)$.  The index $m$ is $N_x$ if the system
includes the Dirichlet condition as a separate equation and $N_x-1$
otherwise.  The matrix $A(u^{-})$ is tridiagonal, so the solution
procedure is to fill a tridiagonal matrix data structure and the
right-hand side vector with the right numbers and call a Gaussian
elimination routine for tridiagonal linear systems.

### Mesh with two cells

mathcal{I}_t helps on the understanding of the details to write out all the
mathematics in a specific
case with a small mesh, say just two cells ($N_x=2$). We use $u^{-}_i$
for the $i$-th component in $u^{-}$.

The starting point is the basic expressions for the
nonlinear equations at mesh point $i=0$ and $i=1$:

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:2x2:x0"></div>

$$
\begin{equation}
A_{0,-1}u_{-1} + A_{0,0}u_0 + A_{0,1}u_1 = b_0,
\label{nonlin:alglevel:1D:fd:2x2:x0} \tag{30}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:2x2:x1"></div>

$$
\begin{equation}  
A_{1,0}u_{0} + A_{1,1}u_1 + A_{1,2}u_2 = b_1\thinspace .
\label{nonlin:alglevel:1D:fd:2x2:x1} \tag{31}
\end{equation}
$$

Equation ([30](#nonlin:alglevel:1D:fd:2x2:x0)) written out reads

$$
\begin{align*}
\frac{1}{2\Delta x^2}(& -(\dfc(u_{-1}) + \dfc(u_{0}))u_{-1}\, +\\ 
& (\dfc(u_{-1}) + 2\dfc(u_{0}) + \dfc(u_{1}))u_0\, -\\ 
& (\dfc(u_{0}) + \dfc(u_{1})))u_1 + au_0
=f(u_0)\thinspace .
\end{align*}
$$

We must then replace $u_{-1}$ by
([29](#nonlin:alglevel:1D:fd:Neumann:x0:um1)).
With Picard iteration we get
<!-- u_{-1} = u_1 -\frac{2\Delta x}{\dfc(u_0)}C -->

$$
\begin{align*}
\frac{1}{2\Delta x^2}(& -(\dfc(u^-_{-1}) + 2\dfc(u^-_{0})
+ \dfc(u^-_{1}))u_1\, +\\ 
&(\dfc(u^-_{-1}) + 2\dfc(u^-_{0}) + \dfc(u^-_{1}))u_0
 + au_0\\ 
&=f(u^-_0) -
\frac{1}{\dfc(u^-_0)\Delta x}(\dfc(u^-_{-1}) + \dfc(u^-_{0}))C,
\end{align*}
$$

where

$$
u^-_{-1} = u_1^- -\frac{2\Delta x}{\dfc(u^-_0)}C\thinspace .
$$

Equation ([31](#nonlin:alglevel:1D:fd:2x2:x1)) contains the unknown $u_2$
for which we have a Dirichlet condition. In case we omit the
condition as a separate equation, ([31](#nonlin:alglevel:1D:fd:2x2:x1))
with Picard iteration becomes

$$
\begin{align*}
\frac{1}{2\Delta x^2}(&-(\dfc(u^-_{0}) + \dfc(u^-_{1}))u_{0}\, + \\ 
&(\dfc(u^-_{0}) + 2\dfc(u^-_{1}) + \dfc(u^-_{2}))u_1\, -\\ 
&(\dfc(u^-_{1}) + \dfc(u^-_{2})))u_2 + au_1
=f(u^-_1)\thinspace .
\end{align*}
$$

We must now move the $u_2$ term to the right-hand side and replace all
occurrences of $u_2$ by $D$:

$$
\begin{align*}
\frac{1}{2\Delta x^2}(&-(\dfc(u^-_{0}) + \dfc(u^-_{1}))u_{0}\, +\\ 
& (\dfc(u^-_{0}) + 2\dfc(u^-_{1}) + \dfc(D)))u_1 + au_1\\ 
&=f(u^-_1) + \frac{1}{2\Delta x^2}(\dfc(u^-_{1}) + \dfc(D))D\thinspace .
\end{align*}
$$

The two equations can be written as a $2\times 2$ system:

$$
\left(\begin{array}{cc}
B_{0,0}& B_{0,1}\\ 
B_{1,0} & B_{1,1}
\end{array}\right)
\left(\begin{array}{c}
u_0\\ 
u_1
\end{array}\right)
=
\left(\begin{array}{c}
d_0\\ 
d_1
\end{array}\right),
$$

where

<!-- Equation labels as ordinary links -->
<div id="_auto10"></div>

$$
\begin{equation}
B_{0,0} =\frac{1}{2\Delta x^2}(\dfc(u^-_{-1}) + 2\dfc(u^-_{0}) + \dfc(u^-_{1}))
+ a,
\label{_auto10} \tag{32}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto11"></div>

$$
\begin{equation}  
B_{0,1} =
-\frac{1}{2\Delta x^2}(\dfc(u^-_{-1}) + 2\dfc(u^-_{0})
+ \dfc(u^-_{1})),
\label{_auto11} \tag{33}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto12"></div>

$$
\begin{equation}  
B_{1,0} =
-\frac{1}{2\Delta x^2}(\dfc(u^-_{0}) + \dfc(u^-_{1})),
\label{_auto12} \tag{34}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto13"></div>

$$
\begin{equation}  
B_{1,1} =
\frac{1}{2\Delta x^2}(\dfc(u^-_{0}) + 2\dfc(u^-_{1}) + \dfc(D)) + a,
\label{_auto13} \tag{35}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto14"></div>

$$
\begin{equation}  
d_0 =
f(u^-_0) -
\frac{1}{\dfc(u^-_0)\Delta x}(\dfc(u^-_{-1}) + \dfc(u^-_{0}))C,
\label{_auto14} \tag{36}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto15"></div>

$$
\begin{equation}  
d_1 = f(u^-_1) + \frac{1}{2\Delta x^2}(\dfc(u^-_{1}) + \dfc(D))D\thinspace .
\label{_auto15} \tag{37}
\end{equation}
$$

The system with the Dirichlet condition becomes

$$
\left(\begin{array}{ccc}
B_{0,0}& B_{0,1} & 0\\ 
B_{1,0} & B_{1,1} & B_{1,2}\\ 
0  & 0 & 1
\end{array}\right)
\left(\begin{array}{c}
u_0\\ 
u_1\\ 
u_2
\end{array}\right)
=
\left(\begin{array}{c}
d_0\\ 
d_1\\ 
D
\end{array}\right),
$$

with

<!-- Equation labels as ordinary links -->
<div id="_auto16"></div>

$$
\begin{equation}
B_{1,1} =
\frac{1}{2\Delta x^2}(\dfc(u^-_{0}) + 2\dfc(u^-_{1}) + \dfc(u_2)) + a,
\label{_auto16} \tag{38}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto17"></div>

$$
\begin{equation}  
B_{1,2} = -
\frac{1}{2\Delta x^2}(\dfc(u^-_{1}) + \dfc(u_2))),
\label{_auto17} \tag{39}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<div id="_auto18"></div>

$$
\begin{equation}  
d_1 = f(u^-_1)\thinspace .
\label{_auto18} \tag{40}
\end{equation}
$$

Other entries are as in the $2\times 2$ system.


### Newton's method

The Jacobian must be derived in order to use Newton's method. Here it means
that we need to differentiate $F(u)=A(u)u - b(u)$ with respect to
the unknown parameters
$u_0,u_1,\ldots,u_m$ ($m=N_x$ or $m=N_x-1$, depending on whether the
Dirichlet condition is included in the nonlinear system $F(u)=0$ or not).
Nonlinear equation number $i$ has the structure

$$
F_i = A_{i,i-1}(u_{i-1},u_i)u_{i-1} +
A_{i,i}(u_{i-1},u_i,u_{i+1})u_i +
A_{i,i+1}(u_i, u_{i+1})u_{i+1} - b_i(u_i)\thinspace .
$$

Computing the Jacobian requires careful differentiation. For example,

$$
\begin{align*}
\frac{\partial}{\partial u_i}(A_{i,i}(u_{i-1},u_i,u_{i+1})u_i) &=
\frac{\partial A_{i,i}}{\partial u_i}u_i + A_{i,i}
\frac{\partial u_i}{\partial u_i}\\ 
&=
\frac{\partial}{\partial u_i}(
\frac{1}{2\Delta x^2}(\dfc(u_{i-1}) + 2\dfc(u_{i})
+\dfc(u_{i+1})) + a)u_i +\\ 
&\quad\frac{1}{2\Delta x^2}(\dfc(u_{i-1}) + 2\dfc(u_{i})
+\dfc(u_{i+1})) + a\\ 
&= \frac{1}{2\Delta x^2}(2\dfc^\prime (u_i)u_i
+\dfc(u_{i-1}) + 2\dfc(u_{i})
+\dfc(u_{i+1})) + a\thinspace .
\end{align*}
$$

The complete Jacobian becomes

$$
\begin{align*}
J_{i,i} &= \frac{\partial F_i}{\partial u_i}
= \frac{\partial A_{i,i-1}}{\partial u_i}u_{i-1}
+ \frac{\partial A_{i,i}}{\partial u_i}u_i
+ A_{i,i}
+ \frac{\partial A_{i,i+1}}{\partial u_i}u_{i+1}
- \frac{\partial b_i}{\partial u_{i}}\\ 
&=
\frac{1}{2\Delta x^2}(
-\dfc^{\prime}(u_i)u_{i-1}
+2\dfc^{\prime}(u_i)u_{i}
+\dfc(u_{i-1}) + 2\dfc(u_i) + \dfc(u_{i+1})) +\\ 
&\quad a
-\frac{1}{2\Delta x^2}\dfc^{\prime}(u_{i})u_{i+1}
- b^{\prime}(u_i),\\ 
J_{i,i-1} &= \frac{\partial F_i}{\partial u_{i-1}}
= \frac{\partial A_{i,i-1}}{\partial u_{i-1}}u_{i-1}
+ A_{i-1,i}
+ \frac{\partial A_{i,i}}{\partial u_{i-1}}u_i
- \frac{\partial b_i}{\partial u_{i-1}}\\ 
&=
\frac{1}{2\Delta x^2}(
-\dfc^{\prime}(u_{i-1})u_{i-1} - (\dfc(u_{i-1}) + \dfc(u_i))
+ \dfc^{\prime}(u_{i-1})u_i),\\ 
J_{i,i+1} &= \frac{\partial A_{i,i+1}}{\partial u_{i-1}}u_{i+1}
+ A_{i+1,i} +
\frac{\partial A_{i,i}}{\partial u_{i+1}}u_i
- \frac{\partial b_i}{\partial u_{i+1}}\\ 
&=\frac{1}{2\Delta x^2}(
-\dfc^{\prime}(u_{i+1})u_{i+1} - (\dfc(u_{i}) + \dfc(u_{i+1}))
+ \dfc^{\prime}(u_{i+1})u_i)
\thinspace .
\end{align*}
$$

The explicit expression for nonlinear equation number $i$,
$F_i(u_0,u_1,\ldots)$, arises from moving the $f(u_i)$ term in
([27](#nonlin:alglevel:1D:fd:deq)) to the left-hand side:

$$
F_i = -\frac{1}{2\Delta x^2}
\left((\dfc(u_i)+\dfc(u_{i+1}))(u_{i+1}-u_i) -
(\dfc(u_{i-1})+\dfc(u_{i}))(u_{i}-u_{i-1})\right)\nonumber
$$

<!-- Equation labels as ordinary links -->
<div id="nonlin:alglevel:1D:fd:deq2"></div>

$$
\begin{equation}  
\qquad\qquad + au_i - f(u_i) = 0\thinspace .
\label{nonlin:alglevel:1D:fd:deq2} \tag{41}
\end{equation}
$$

At the boundary point $i=0$, $u_{-1}$ must be replaced using
the formula ([29](#nonlin:alglevel:1D:fd:Neumann:x0:um1)).
When the Dirichlet condition at $i=N_x$ is not a part of the
equation system, the last equation $F_m=0$ for $m=N_x-1$
involves the quantity $u_{N_x-1}$ which must be replaced by $D$.
If $u_{N_x}$ is treated as an unknown in the system, the
last equation $F_m=0$ has $m=N_x$ and reads

$$
F_{N_x}(u_0,\ldots,u_{N_x}) = u_{N_x} - D = 0\thinspace .
$$

Similar replacement of $u_{-1}$ and $u_{N_x}$ must be done in
the Jacobian for the first and last row. When $u_{N_x}$
is included as an unknown, the last row in the Jacobian
must help implement the condition $\delta u_{N_x}=0$, since
we assume that $u$ contains the right Dirichlet value
at the beginning of the iteration ($u_{N_x}=D$), and then
the Newton update should be zero for $i=0$, i.e., $\delta u_{N_x}=0$.
This also forces the right-hand side to be $b_i=0$, $i=N_x$.

We have seen, and can see from the present example, that the
linear system in Newton's method contains all the terms present
in the system that arises in the Picard iteration method.
The extra terms in Newton's method can be multiplied by a factor
such that it is easy to program one linear system and set this
factor to 0 or 1 to generate the Picard or Newton system.

<!-- Remark: Neumann cond at x=L and Dirichlet at x=0 leads to different -->
<!-- numbering of unknowns and u at mesh points. Must address this -->
<!-- in a remark and treat it properly in diffu. -->

f08d1bae967493b89be4797f48b44518194ba211
ipynb
Jupyter Notebook
fdm-devito-notebooks/05_nonlin/nonlin_pde1D.ipynb
devitocodes/devito_book
30405c3d440a1f89df69594fd0704f69650c1ded
CC-BY-4.0
2020-07-17T13:19:15.000Z
2021-03-27T05:21:09.000Z
fdm-devito-notebooks/05_nonlin/nonlin_pde1D.ipynb
devitocodes/devito_book
30405c3d440a1f89df69594fd0704f69650c1ded
CC-BY-4.0
2020-07-14T15:38:52.000Z
2020-09-25T11:54:59.000Z
fdm-devito-notebooks/05_nonlin/nonlin_pde1D.ipynb
devitocodes/devito_book
30405c3d440a1f89df69594fd0704f69650c1ded
CC-BY-4.0
2021-03-27T05:21:14.000Z
2021-03-27T05:21:14.000Z
Qwen/Qwen-72B
1. YES
2. YES

