Discretization of 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
The problem (47) arises from the stationary limit of a diffusion equation,
as \(t\rightarrow\infty\) and \(\partial u/\partial t\rightarrow 0\). Alternatively, the problem (47) arises at each time level from implicit time discretization of (48). For example, a Backward Euler scheme for (48) leads to
Introducing \(u(x)\) for \(u^n(x)\), \(u^{(1)}\) for \(u^{n-1}\), and defining \(f(u)\) in (47) to be \(f(u)\) in (49) plus \(u^{n-1}/\Delta t\), gives (47) with \(a=1/\Delta t\).
Finite difference discretization¶
Since the technical steps in finite difference discretization in space are so much simpler than the steps in the finite element method, we start with finite difference to illustrate the concept of handling this nonlinear problem and minimize the spatial discretization details.
The nonlinearity in the differential equation (47) 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:
Writing this out for a uniform mesh with points \(x_i=i\Delta x\), \(i=0,\ldots,N_x\), leads to
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=D\). 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 on \(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:
Equation (50) with the latter approximation then looks like
or written more compactly,
At mesh point \(i=0\) we have the boundary condition \({\alpha}(u)u^{\prime}=C\), which is discretized by
meaning
The fictitious value \(u_{-1}\) can be eliminated with the aid of (51) for \(i=0\). Formally, (51) should be solved with respect to \(u_{i-1}\) and that value (for \(i=0\)) should be inserted in (52), 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 (52) after every Picard or Newton iteration. Such an approach means that we use a known \(u_{-1}\) value in (51) from the previous iteration.
Solution of algebraic equations¶
The structure of the equation system¶
The nonlinear algebraic equations (51) are of the form \(A(u)u = b(u)\) with
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 (52):
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
Equation (54) written out reads
We must then replace \(u_{-1}\) by (53). With Picard iteration we get
where
Equation (55) contains the unknown \(u_2\) for which we have a Dirichlet condition. In case we omit the condition as a separate equation, (55) with Picard iteration becomes
We must now move the \(u_2\) term to the right-hand side and replace all occurrences of \(u_2\) by \(D\):
The two equations can be written as a \(2\times 2\) system:
where
The system with the Dirichlet condition becomes
with
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
Computing the Jacobian requires careful differentiation. For example,
The complete Jacobian becomes
The explicit expression for nonlinear equation number \(i\), \(F_i(u_0,u_1,\ldots)\), arises from moving the \(f(u_i)\) term in (51) to the left-hand side:
At the boundary point \(i=0\), \(u_{-1}\) must be replaced using the formula (53). 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
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.
Galerkin-type discretization¶
For a Galerkin-type discretization, which may be developed into a finite element method, we first need to derive the variational problem. Let \(V\) be an appropriate function space with basis functions \(\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}\). Because of the Dirichlet condition at \(x=L\) we require \({\psi}_i(L)=0\), \(i\in{\mathcal{I}_s}\). The approximate solution is written as \(u = D + \sum_{j\in{\mathcal{I}_s}}c_j{\psi}_j\), where the term \(D\) can be viewed as a boundary function needed to implement the Dirichlet condition \(u(L)=D\).
Using Galerkin’s method, we multiply the differential equation by any \(v\in V\) and integrate terms with second-order derivatives by parts:
The Neumann condition at the boundary \(x=0\) is inserted in the boundary term:
(Recall that since \({\psi}_i(L)=0\), any linear combination \(v\) of the basis functions also vanishes at \(x=L\): \(v(L)=0\).) The variational problem is then: find \(u\in V\) such that
To derive the algebraic equations, we note that \(\forall v\in V\) is equivalent with \(v={\psi}_i\) for \(i\in{\mathcal{I}_s}\). Setting \(u=D+\sum_jc_j{\psi}_j\) and sorting terms results in the linear system
Fundamental integration problem¶
Methods that use the Galerkin or weighted residual method face a fundamental difficulty in nonlinear problems: how can we integrate a terms like \(\int_0^L {\alpha}(\sum_{k}c_k{\psi}_k){\psi}_i^{\prime}{\psi}_j^{\prime}{\, \mathrm{d}x}\) and \(\int_0^L f(\sum_{k}c_k{\psi}_k){\psi}_i{\, \mathrm{d}x}\) when we do not know the \(c_k\) coefficients in the argument of the \({\alpha}\) and \(f\) functions? We can resort to numerical integration, provided an approximate \(\sum_kc_k{\psi}_k\) can be used for the argument \(u\) in \(f\) and \({\alpha}\). This is the approach used in computer programs.
However, if we want to look more mathematically into the structure of the algebraic equations generated by the finite element method in nonlinear problems, and compare such equations with those arising in the finite difference method, we need techniques that enable integration of expressions like \(\int_0^L f(\sum_{k}c_k{\psi}_k){\psi}_i{\, \mathrm{d}x}\) by hand. Two such techniques will be shown: the group finite element and numerical integration based on the nodes only. Both techniques are approximate, but they allow us to see the difference equations in the finite element method. The details are worked out in Appendix Appendix: Symbolic nonlinear finite element equations. Some readers will prefer to dive into these symbolic calculations to gain more understanding of nonlinear finite element equations, while others will prefer to continue with computational algorithms (in the next two sections) rather than analysis.
Picard iteration defined from the variational form¶
Consider the problem (47) with the corresponding variational form (66). Our aim is to define a Picard iteration based on this variational form without any attempt to compute integrals symbolically as in the previous three sections. The idea in Picard iteration is to use a previously computed \(u\) value in the nonlinear functions \({\alpha}(u)\) and \(f(u)\). Let \(u^{-}\) be the available approximation to \(u\) from the previous Picard iteration. The linearized variational form for Picard iteration is then
This is a linear problem \(a(u,v)=L(v)\) with bilinear and linear forms
Make sure to distinguish the coefficient \(a\) in \(auv\) from the differential equation from the \(a\) in the abstract bilinear form notation \(a(\cdot,\cdot)\).
The linear system associated with (68) is computed in the standard way. Technically, we are back to solving \(-({\alpha}(x)u^{\prime})^{\prime} + au=f(x)\). The unknown \(u\) is sought on the form \(u = B(x) + \sum_{j\in{\mathcal{I}_s}}c_j{\psi}_j\), with \(B(x)=D\) and \({\psi}_i = {\varphi}_{\nu(i)}\), \(\nu(i)=i\), and \({\mathcal{I}_s} = \{0,1,\ldots,N=N_n-2\}\).
Newton’s method defined from the variational form¶
Application of Newton’s method to the nonlinear variational form (66) arising from the problem (47) requires identification of the nonlinear algebraic equations \(F_i=0\). Although we originally denoted the unknowns in nonlinear algebraic equations by \(u_0,\ldots,u_N\), it is in the present context most natural to have the unknowns as \(c_0,\ldots,c_N\) and write
and define the Jacobian as \(J_{i,j}=\partial F_i/\partial c_j\) for \(i,j\in{\mathcal{I}_s}\).
The specific form of the equations \(F_i=0\) follows from the variational form
by choosing \(v={\psi}_i\), \(i\in{\mathcal{I}_s}\), and setting \(u=\sum_{j\in{\mathcal{I}_s}}c_j{\psi}_j\), maybe with a boundary function to incorporate Dirichlet conditions.
With \(v={\psi}_i\) we get
In the differentiations leading to the Jacobian we will frequently use the results
The derivation of the Jacobian of (69) goes as
One must be careful about the prime symbol as differentiation
In \({\alpha}^{\prime}\) the derivative is with respect to the independent variable in the \({\alpha}\) function, and that is \(u\), so
while in \(u^{\prime}\) the differentiation is with respect to \(x\), so
Similarly, \(f\) is a function of \(u\), so \(f^{\prime}\) means \(df/du\).
When calculating the right-hand side vector \(F_i\) and the coefficient matrix \(J_{i,j}\) in the linear system to be solved in each Newton iteration, one must use a previously computed \(u\), denoted by \(u^{-}\), for the symbol \(u\) in (69) and (70). With this notation we have
These expressions can be used for any basis \(\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}\). Choosing finite element functions for \({\psi}_i\), one will normally want to compute the integral contribution cell by cell, working in a reference cell. To this end, we restrict the integration to one cell and transform the cell to \([-1,1]\). The most recently computed approximation \(u^{-}\) to \(u\) becomes \(\tilde u^{-}=\sum_t\tilde u^{-1}_t{\tilde{\varphi}}_t(X)\) over the reference element, where \(\tilde u^{-1}_t\) is the value of \(u^{-}\) at global node (or degree of freedom) \(q(e,t)\) corresponding to the local node \(t\) (or degree of freedom). The formulas (71) and (72) then change to
with \(r,s\in{I_d}\) runs over the local degrees of freedom.
Many finite element programs require the user to provide \(F_i\) and \(J_{i,j}\). Some programs, like FEniCS, are capable of automatically deriving \(J_{i,j}\) if \(F_i\) is specified.
Dirichlet conditions¶
Incorporation of the Dirichlet values by assembling contributions from all degrees of freedom and then modifying the linear system can obviously be applied to Picard iteration as that method involves a standard linear system. In the Newton system, however, the unknown is a correction \(\delta u\) to the solution. Dirichlet conditions are implemented by inserting them in the initial guess \(u^{-}\) for the Newton iteration and implementing \(\delta u_i =0\) for all known degrees of freedom. The manipulation of the linear system follows exactly the algorithm in the linear problems, the only difference being that the known values are zero.
