1D stationary nonlinear differential equations

The section Linearization at the differential equation level 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 Systems of nonlinear algebraic equations. This latter approach will now be described in detail.

We shall work with the 1D problem

$$\tag{569} -({\alpha}(u)u^{\prime})^{\prime} + au = f(u),\quad x\in (0,L), \quad {\alpha}(u(0))u^{\prime}(0) = C,\ u(L)=D {\thinspace .}$$

The problem (569) arises from the stationary limit of a diffusion equation,

$$\tag{570} \frac{\partial u}{\partial t} = \frac{\partial}{\partial x}\left( \alpha(u)\frac{\partial u}{\partial x}\right) - au + f(u),$$

as $t\rightarrow\infty$ and $\partial u/\partial t\rightarrow 0$. Alternatively, the problem (569) arises at each time level from implicit time discretization of (570). For example, a Backward Euler scheme for (570) leads to

$$\tag{571} \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 .}$$

Introducing $u(x)$ for $u^n(x)$, $u^{(1)}$ for $u^{n-1}$, and defining $f(u)$ in (569) to be $f(u)$ in (571) plus $u^{n-1}/\Delta t$, gives (569) with $a=1/\Delta t$.

Finite difference discretization

The nonlinearity in the differential equation (569) poses no more difficulty than a variable coefficient, as in the term $({\alpha}(x)u^{\prime})^{\prime}$. We can therefore use a standard finite difference approach to discretizing the Laplace term with a variable coefficient:

$$[-D_x{\alpha} 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

$$\tag{572} -\frac{1}{\Delta x^2} \left({\alpha}_{i+\frac{1}{2}}(u_{i+1}-u_i) - {\alpha}_{i-\frac{1}{2}}(u_{i}-u_{i-1})\right) + au_i = f(u_i){\thinspace .}$$

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 $({\alpha}(x)u^{\prime})^{\prime}$ and $f(x)$ is that now ${\alpha}$ and $f$ are functions of $u$ and not only of $x$: $({\alpha}(u(x))u^{\prime})^{\prime}$ and $f(u(x))$.

The quantity ${\alpha}_{i+\frac{1}{2}}$, evaluated between two mesh points, needs a comment. Since ${\alpha}$ depends on $u$ and $u$ is only known at the mesh points, we need to express ${\alpha}_{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 ${\alpha}$ features large jumps. There are two choices of arithmetic means:

$$\tag{573} {\alpha}_{i+\frac{1}{2}} \approx {\alpha}(\frac{1}{2}(u_i + u_{i+1}) = [{\alpha}(\overline{u}^x)]^{i+\frac{1}{2}},$$

$$\tag{574} {\alpha}_{i+\frac{1}{2}} \approx \frac{1}{2}({\alpha}(u_i) + {\alpha}(u_{i+1})) = [\overline{{\alpha}(u)}^x]^{i+\frac{1}{2}}$$

Equation (572) with the latter approximation then looks like

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

$$\tag{575} \qquad\qquad + au_i = f(u_i),$$

or written more compactly,

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

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

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

meaning

$$\tag{576} {\alpha}(u_0)\frac{u_{1} - u_{-1}}{2\Delta x} = C{\thinspace .}$$

The fictitious value $u_{-1}$ can be eliminated with the aid of (575) for $i=0$. Formally, (575) should be solved with respect to $u_{i-1}$ and that value (for $i=0$) should be inserted in (576), 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 (576) after every Picard or Newton iteration. Such an approach means that we use a known $u_{-1}$ value in (575) from the previous iteration.

Solution of algebraic equations

The structure of the equation system

The nonlinear algebraic equations (575) are of the form $A(u)u = b(u)$ with

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

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 (576):

$$\tag{577} u_{-1} = u_1 -\frac{2\Delta x}{{\alpha}(u_0)}C{\thinspace .}$$

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$. It 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 Systems of nonlinear algebraic equations. 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

It 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$ are

$$\tag{578} A_{0,-1}u_{-1} + A_{0,0}u_0 + A_{0,1}u_1 = b_0,$$

$$\tag{579} A_{1,0}u_{0} + A_{1,1}u_1 + A_{1,2}u_2 = b_1{\thinspace .}$$

Equation (578) written out reads

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

We must then replace $u_{-1}$ by (577). With Picard iteration we get

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

where

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

Equation (579) contains the unknown $u_2$ for which we have a Dirichlet condition. In case we omit the condition as a separate equation, (579) with Picard iteration becomes

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

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

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

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

$$\begin{split}\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),\end{split}$$

where

$$\tag{580} B_{0,0} =\frac{1}{2\Delta x^2}({\alpha}(u^-_{-1}) + 2{\alpha}(u^-_{0}) + {\alpha}(u^-_{1})) + a$$

$$\tag{581} B_{0,1} = -\frac{1}{2\Delta x^2}({\alpha}(u^-_{-1}) + 2{\alpha}(u^-_{0}) + {\alpha}(u^-_{1})),$$

$$\tag{582} B_{1,0} = -\frac{1}{2\Delta x^2}({\alpha}(u^-_{0}) + {\alpha}(u^-_{1})),$$

$$\tag{583} B_{1,1} = \frac{1}{2\Delta x^2}({\alpha}(u^-_{0}) + 2{\alpha}(u^-_{1}) + {\alpha}(D)) + a,$$

$$\tag{584} d_0 = f(u^-_0) - \frac{1}{{\alpha}(u^-_0)\Delta x}({\alpha}(u^-_{-1}) + {\alpha}(u^-_{0}))C,$$

$$\tag{585} d_1 = f(u^-_1) + \frac{1}{2\Delta x^2}({\alpha}(u^-_{1}) + {\alpha}(D))D{\thinspace .}$$

The system with the Dirichlet condition becomes

$$\begin{split}\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),\end{split}$$

with

$$\tag{586} B_{1,1} = \frac{1}{2\Delta x^2}({\alpha}(u^-_{0}) + 2{\alpha}(u^-_{1}) + {\alpha}(u_2)) + a,$$

$$\tag{587} B_{1,2} = - \frac{1}{2\Delta x^2}({\alpha}(u^-_{1}) + {\alpha}(u_2))),$$

$$\tag{588} d_1 = f(u^-_1){\thinspace .}$$

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{split}\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}({\alpha}(u_{i-1}) + 2{\alpha}(u_{i}) +{\alpha}(u_{i+1})) + a)u_i +\\ &\quad\frac{1}{2\Delta x^2}({\alpha}(u_{i-1}) + 2{\alpha}(u_{i}) +{\alpha}(u_{i+1})) + a\\ &= \frac{1}{2\Delta x^2}(2{\alpha}^\prime (u_i)u_i +{\alpha}(u_{i-1}) + 2{\alpha}(u_{i}) +{\alpha}(u_{i+1})) + a{\thinspace .} \end{align*}\end{split}$$

The complete Jacobian becomes

$$\begin{split}\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}( -{\alpha}^{\prime}(u_i)u_{i-1} +2{\alpha}^{\prime}(u_i)u_{i} +{\alpha}(u_{i-1}) + 2{\alpha}(u_i) + {\alpha}(u_{i+1})) +\\ &\quad a -\frac{1}{2\Delta x^2}{\alpha}^{\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}( -{\alpha}^{\prime}(u_{i-1})u_{i-1} - ({\alpha}(u_{i-1}) + {\alpha}(u_i)) + {\alpha}^{\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}( -{\alpha}^{\prime}(u_{i+1})u_{i+1} - ({\alpha}(u_{i}) + {\alpha}(u_{i+1})) + {\alpha}^{\prime}(u_{i+1})u_i) {\thinspace .} \end{align*}\end{split}$$

The explicit expression for nonlinear equation number $i$, $F_i(u_0,u_1,\ldots)$, arises from moving the $f(u_i)$ term in (575) to the left-hand side:

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

$$\tag{589} \qquad\qquad + au_i - f(u_i) = 0{\thinspace .}$$

At the boundary point $i=0$, $u_{-1}$ must be replaced using the formula (577). 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.