The integrals in nonlinear finite element equations are computed by numerical integration rules in computer programs, so the formulas for the variational form are directly transferred to numbers. It is also of interest to understand the nature of the system of difference equations that arises from the finite element method in nonlinear problems and to compare with corresponding expressions arising from finite difference discretization. We shall dive into this problem here. To see the structure of the difference equations implied by the finite element method, we have to find symbolic expressions for the integrals, and this is extremely difficult since the integrals involve the unknown function in nonlinear problems. However, there are some techniques that allow us to approximate the integrals and work out symbolic formulas that can be compared with their finite difference counterparts.
We shall address the 1D model problem (50) from the beginning of the section Discretization of 1D stationary nonlinear differential equations. The finite difference discretization is shown in the section Finite difference discretization, while the variational form based on Galerkin's method is developed in the section Galerkin-type discretization. We build directly on formulas developed in the latter section.
Introduction of finite element basis functions \( \basphi_i \) means setting $$ \baspsi_i = \basphi_{\nu(i)},\quad i\in\If,$$ where degree of freedom number \( \nu(i) \) in the mesh corresponds to unknown number \( i \) (\( c_i \)). In the present example, we use all the basis functions except the last at \( i=N_n-1 \), i.e., \( \If = \{0,\ldots,N_n-2\} \), and \( \nu(j)=j \). The expansion of \( u \) can be taken as $$ u = D + \sum_{j\in\If} c_j\basphi_{\nu(j)},$$ but it is more common in a finite element context to use a boundary function \( B=\sum_{j\in\Ifb}U_j\basphi_j \), where \( U_j \) is prescribed Dirichlet condition for degree of freedom number \( j \) and \( U_j \) is the corresponding value. $$ u = D\basphi_{N_n-1} + \sum_{j\in\If} c_j\basphi_{\nu(j)}\tp $$ In the general case with \( u \) prescribed as \( U_j \) at some nodes \( j\in\Ifb \), we set $$ u = \sum_{j\in\Ifb} U_j\basphi_j + \sum_{j\in\If}c_j\basphi_{\nu(j)},$$ where \( c_j = u(x_{\nu(j)}) \). That is, \( \nu(j) \) maps unknown number \( j \) to the corresponding node number \( \nu(j) \) such that \( c_j = u(x_{\nu(j)}) \).
Since we already expand \( u \) as \( \sum_j\basphi_j u(\xno{j}) \), we may use the same approximation for other functions as well. For example, $$ f(u)\approx \sum_{j} f(\xno{j})\basphi_j, $$ where \( f(\xno{j}) \) is the value of \( f \) at node \( j \). Since \( f \) is a function of \( u \), \( f(\xno{j})=f(u(\xno{j})) \). Introducing \( u_j \) as a short form for \( u(\xno{j}) \), we can write $$ f(u)\approx \sum_{j} f(u_{j})\basphi_j\tp $$ This approximation is known as the group finite element method or the product approximation technique. The index \( j \) runs over all node numbers in the mesh.
The principal advantages of the group finite element method are two-fold:
Note, however, that an expression like \( \int f(u)\basphi_i\dx \) causes no problems in a computer program as the integral is calculated by numerical integration using an existing approximation of \( u \) in \( f(u) \) such that the integrand can be sampled at any spatial point.
Our aim now is to derive symbolic expressions for the difference equations arising from the finite element method in nonlinear problems and compare the expressions with those arising in the finite difference method. To this end, let us simplify the model problem and set \( a=0 \), \( \dfc=1 \), \( f(u)=u^2 \), and have Neumann conditions at both ends such that we get a very simple nonlinear problem \( -u^{\prime\prime}=u^2 \), \( u'(0)=1 \), \( u'(L)=0 \). The variational form is then $$ \int_0^L u^{\prime}v^{\prime}\dx = \int_0^L u^2v\dx - v(0),\quad\forall v\in V\tp$$ The term with \( u^{\prime}v^{\prime} \) is well known so the only new feature is the term \( \int u^2v\dx \).
To make the distance from finite element equations to finite difference equations as short as possible, we shall substitute \( c_j \) in the sum \( u=\sum_jc_j\basphi_j \) by \( u_j=u(\xno{j}) \) since \( c_j \) is the value of \( u \) at node \( j \). (In the more general case with Dirichlet conditions as well, we have a sum \( \sum_jc_j\basphi_{\nu(j)} \) where \( c_j \) is replaced by \( u(\xno{\nu(j)}) \)). We can then introduce some other counter \( k \) such that it is meaningful to write \( u=\sum_k u_k\basphi_k \), where \( k \) runs over appropriate node numbers.) The quantity \( u_j \) in \( \sum_ju_j\basphi_j \) is the same as \( u \) at mesh point number \( j \) in the finite difference method, which is commonly denoted \( u_j \).
Consider the term \( \int u^2v\dx \) in the variational formulation with \( v=\basphi_i \) and \( u=\sum_k\basphi_ku_k \): $$ \int_0^L (\sum_ku_k\basphi_k)^2\basphi_i\dx\tp$$ Evaluating this integral for P1 elements (see Problem 11: Integrate functions of finite element expansions) results in the expression $$ \frac{h}{12}(u_{i-1}^2 + 2u_i(u_{i-1} + u_{i+1}) + 6u_i^2 + u_{i+1}^2),$$ to be compared with the simple value \( u_i^2 \) that would arise in a finite difference discretization when \( u^2 \) is sampled at mesh point \( x_i \). More complicated \( f(u) \) functions in the integral \( \int_0^L f(u)\basphi_i\dx \) give rise to much more lengthy expressions, if it is possible to carry out the integral symbolically at all.
Let use the group finite element method to derive the terms in the difference equation corresponding to \( f(u) \) in the differential equation. We have $$ \int_0^L f(u)\basphi_i\dx \approx \int_0^L (\sum_j \basphi_jf(u_j))\basphi_i\dx = \sum_j \left(\int_0^L \basphi_i\basphi_j\dx\right) f(u_j)\tp$$ We recognize this expression as the mass matrix \( M \), arising from \( \int\basphi_i\basphi_j\dx \), times the vector \( f=(f(u_0),f(u_1),\ldots,) \): \( Mf \). The associated terms in the difference equations are, for P1 elements, $$ \frac{h}{6}(f(u_{i-1}) + 4f(u_i) + f(u_{i+1}))\tp$$ Occasionally, we want to interpret this expression in terms of finite differences, and to this end a rewrite of this expression is convenient: $$ \frac{h}{6}(f(u_{i-1}) + 4f(u_i) + f(u_{i+1})) = h[f(u) - \frac{h^2}{6}D_xD_x f(u)]_i\tp$$ That is, the finite element treatment of \( f(u) \) (when using a group finite element method) gives the same term as in a finite difference approach, \( f(u_i) \), minus a diffusion term which is the 2nd-order discretization of \( \frac{1}{6}h^2f''(x_i) \).
We may lump the mass matrix through integration with the Trapezoidal rule so that \( M \) becomes diagonal in the finite element method. In that case the \( f(u) \) term in the differential equation gives rise to a single term \( hf(u_i) \), just as in the finite difference method.
Let us reconsider a term \( \int f(u)v\dx \) as treated in the previous section, but now we want to integrate this term numerically. Such an approach can lead to easy-to-interpret formulas if we apply a numerical integration rule that samples the integrand at the node points \( \xno{i} \) only, because at such points, \( \basphi_j(\xno{i})=0 \) if \( j\neq i \), which leads to great simplifications.
The term in question takes the form $$ \int_0^L f(\sum_k u_k\basphi_k)\basphi_i\dx\tp$$ Evaluation of the integrand at a node \( \xno{\ell} \) leads to a collapse of the sum \( \sum_k u_k\basphi_k \) to one term because $$ \sum_k u_k\basphi_k(\xno{\ell}) = u_\ell\tp$$ $$ f(\sum_k u_k\underbrace{\basphi_k(\xno{\ell})}_{\delta_{k\ell}}) \underbrace{\basphi_i(\xno{\ell})}_{\delta_{i\ell}} = f(u_\ell)\delta_{i\ell},$$ where we have used the Kronecker delta: \( \delta_{ij}=0 \) if \( i\neq j \) and \( \delta_{ij}=1 \) if \( i=j \).
Considering the Trapezoidal rule for integration, where the integration points are the nodes, we have $$ \int_0^L f(\sum_k u_k\basphi_k(x))\basphi_i(x)\dx \approx h\sum_{\ell=0}^{N_n-1} f(u_\ell)\delta_{i\ell} - \mathcal{C} = hf(u_i)\tp $$ This is the same representation of the \( f \) term as in the finite difference method. The term \( \mathcal{C} \) contains the evaluations of the integrand at the ends with weight \( \half \), needed to make a true Trapezoidal rule: $$ \mathcal{C} = \frac{h}{2}f(u_0)\basphi_i(0) + \frac{h}{2}f(u_{N_n-1})\basphi_i(L)\tp$$ The answer \( hf(u_i) \) must therefore be multiplied by \( \half \) if \( i=0 \) or \( i=N_n-1 \). Note that \( \mathcal{C}=0 \) for \( i=1,\ldots,N_n-2 \).
One can alternatively use the Trapezoidal rule on the reference cell and assemble the contributions. It is a bit more labor in this context, but working on the reference cell is safer as that approach is guaranteed to handle discontinuous derivatives of finite element functions correctly (not important in this particular example), while the rule above was derived with the assumption that \( f \) is continuous at the integration points.
The conclusion is that it suffices to use the Trapezoidal rule if one wants to derive the difference equations in the finite element method and make them similar to those arising in the finite difference method. The Trapezoidal rule has sufficient accuracy for P1 elements, but for P2 elements one should turn to Simpson's rule.
Turning back to the model problem (50), it remains to calculate the contribution of the \( (\dfc u^{\prime})^{\prime} \) and boundary terms to the difference equations. The integral in the variational form corresponding to \( (\dfc u^{\prime})^{\prime} \) is $$ \int_0^L \dfc(\sum_k c_k\basphi_k)\basphi_i^{\prime}\basphi_j^{\prime}\dx\tp$$ Numerical integration utilizing a value of \( \sum_k c_k\basphi_k \) from a previous iteration must in general be used to compute the integral. Now our aim is to integrate symbolically, as much as we can, to obtain some insight into how the finite element method approximates this term. To be able to derive symbolic expressions, we must either turn to the group finite element method or numerical integration in the node points. Finite element basis functions \( \basphi_i \) are now used.
We set \( \dfc(u)\approx \sum_k\alpha(u_k)\basphi_k \), and then we write $$ \int_0^L \dfc(\sum_k c_k\basphi_k)\basphi_i^{\prime}\basphi_j^{\prime}\dx \approx \sum_k (\underbrace{\int_0^L \basphi_k\basphi_i^{\prime}\basphi_j^{\prime}\dx}_{L_{i,j,k}}) \dfc(u_k) = \sum_k L_{i,j,k}\dfc(u_k)\tp $$ Further calculations are now easiest to carry out in the reference cell. With P1 elements we can compute \( L_{i,j,k} \) for the two \( k \) values that are relevant on the reference cell. Turning to local indices, one gets $$ L_{r,s,t}^{(e)} = \frac{1}{2h}\left(\begin{array}{rr} 1 & -1\\ -1 & 1 \end{array}\right),\quad t=0, 1, $$ where \( r,s,t=0,1 \) are indices over local degrees of freedom in the reference cell (\( i=q(e,r) \), \( j=q(e,s) \), and \( k=q(e,t) \)). The sum \( \sum_k L_{i,j,k}\dfc(u_k) \) at the cell level becomes \( \sum_{t=0}^1 L_{r,s,t}^{(e)}\dfc(\tilde u_t) \), where \( \tilde u_t \) is \( u(\xno{q(e,t)}) \), i.e., the value of \( u \) at local node number \( t \) in cell number \( e \). The element matrix becomes $$ \begin{equation} \half (\dfc(\tilde u_0) + \dfc(\tilde u_{1})) \frac{1}{h}\left(\begin{array}{rr} 1 & -1\\ -1 & 1 \end{array}\right)\tp \tag{101} \end{equation} $$ As usual, we employ a left-to-right numbering of cells and nodes. Row number \( i \) in the global matrix gets contributions from the first row of the element matrix in cell \( i \) and the last row of the element matrix in cell \( i-1 \). In cell number \( i-1 \) the sum \( \dfc(\tilde u_0) + \dfc(\tilde u_{1}) \) corresponds to \( \dfc(u_{i-1}) + \dfc(u_i) \). The same sum becomes \( \dfc(u_{i}) + \dfc(u_{i+1}) \) in cell number \( i \). We can with this insight assemble the contributions to row number \( i \) in the global matrix: $$ \frac{1}{2h}(-(\dfc(u_{i-1}) + \dfc(u_i)),\quad \dfc(u_{i-1}) + 2\dfc(u_i) + \dfc(u_{i+1}),\quad -(\dfc(u_{i}) + \dfc(u_{i+1})))\tp $$ Multiplying by the vector of unknowns \( u_i \) results in a formula that can be arranged to $$ \begin{equation} -\frac{1}{h}(\half(\dfc(u_i) + \dfc(u_{i+1}))(u_{i+1}-u_i) - \half(\dfc(u_{i-1}) + \dfc(u_{i}))(u_{i}-u_{i-1})), \tag{102} \end{equation} $$ which is nothing but the standard finite difference discretization of \( -(\dfc(u)u^{\prime})^{\prime} \) with an arithmetic mean of \( \dfc(u) \) (and the usual factor \( h \) because of the integration in the finite element method).
Instead of using the group finite element method and exact integration we can turn to the Trapezoidal rule for computing \( \int_0^L \dfc(\sum_k u_k\basphi_k)\basphi_i^{\prime}\basphi_j^{\prime}\dx \), again at the cell level since that is most convenient when we deal with discontinuous functions \( \basphi_i' \): $$ \begin{align} \int_{-1}^1 \alpha(\sum_t\tilde u_t\refphi_t)\refphi_r'\refphi_s'\frac{h}{2}dX &= \int_{-1}^1 \dfc(\sum_{t=0}^1 \tilde u_t\refphi_t)\frac{2}{h}\frac{d\refphi_r}{dX} \frac{2}{h}\frac{d\refphi_s}{dX}\frac{h}{2}dX\nonumber\\ & = \frac{1}{2h}(-1)^r(-1)^s \int_{-1}^1 \dfc(\sum_{t=0}^1 u_t\refphi_t(X))dX \nonumber\\ & \approx \frac{1}{2h}(-1)^r(-1)^s\dfc ( \sum_{t=0}^1\refphi_t(-1)\tilde u_t) + \dfc(\sum_{t=0}^1\refphi_t(1)\tilde u_t) \nonumber\\ & = \frac{1}{2h}(-1)^r(-1)^s(\dfc(\tilde u_0) + \dfc(\tilde u_1))\tp \tag{103} \end{align} $$ The element matrix in (103) is identical to the one in (101), showing that the group finite element method and Trapezoidal integration are equivalent with a standard finite discretization of a nonlinear Laplace term \( (\dfc(u)u^{\prime})^{\prime} \) using an arithmetic mean for \( \dfc \): \( [D_x\overline{x}D_xu]_i \).
The term \( \int auv\dx \) in the variational form is linear and gives these terms in the algebraic equations: $$ \frac{ah}{6}(u_{i-1} + 4u_i + u_{i+1}) = ah[u - \frac{h^2}{6}D_xD_x u]_i\tp$$ The final term in the variational form is the Neumann condition at the boundary: \( Cv(0)=C\basphi_i(0) \). With a left-to-right numbering only \( i=0 \) will give a contribution \( Cv(0)=C\delta_{i0} \) (since \( \basphi_i(0)\neq 0 \) only for \( i=0 \)).
As we now have explicit expressions for the nonlinear difference equations also in the finite element method, a Picard or Newton method can be defined as shown for the finite difference method. However, our effort in deriving symbolic forms of the difference equations in the finite element method was motivated by a desire to see how nonlinear terms in differential equations make the finite element and difference method different. For practical calculations in computer programs we apply numerical integration, normally the more accurate Gauss-Legendre quadrature rules, to the integrals directly. This allows us to easily evaluate the nonlinear algebraic equations for a given numerical approximation of \( u \) (here denoted \( u^{-} \)). To solve the nonlinear algebraic equations we need to apply the Picard iteration method or Newton's method to the variational form directly, as shown next.