We have to take special actions to incorporate Dirichlet conditions, such as \( u(L)=D \), into the computational procedures. The present section outlines alternative, yet mathematically equivalent, methods.
In the section Boundary function we introduce a boundary function \( B(x) \) to deal with nonzero Dirichlet boundary conditions for \( u \). The construction of such a function is not always trivial, especially not in multiple dimensions. However, a simple and general construction idea exists when the basis functions have the property
$$ \begin{equation*} \basphi_i(\xno{j}) = \delta_{ij},\quad \delta_{ij} = \left\lbrace\begin{array}{ll} 1, & i=j,\\ 0, & i\neq j, \end{array}\right. \end{equation*} $$ where \( \xno{j} \) is a boundary point. Examples on such functions are the Lagrange interpolating polynomials and finite element functions.
Suppose now that \( u \) has Dirichlet boundary conditions at nodes with numbers \( i\in\Ifb \). For example, \( \Ifb = \{0,N_n\} \) in a 1D mesh with node numbering from left to right. Let \( U_i \) be the corresponding prescribed values of \( u(\xno{i}) \). We can then, in general, use
$$ \begin{equation} B(x) = \sum_{j\in\Ifb} U_j\basphi_j(x)\tp \end{equation} $$ It is easy to verify that \( B(\xno{i})= \sum_{j\in\Ifb} U_j\basphi_j(\xno{i})=U_i \).
The unknown function can then be written as
$$ \begin{equation} u(x) = \sum_{j\in\Ifb} U_j\basphi_j(x) + \sum_{j\in\If}c_j\basphi_{\nu(j)}, \end{equation} $$ where \( \nu(j) \) maps unknown number \( j \) in the equation system to node \( \nu(j) \). We can easily show that with this \( u \), a Dirichlet condition \( u(\xno{k})=U_k \) is fulfilled:
$$ u(\xno{k}) = \sum_{j\in\Ifb} U_j\underbrace{\basphi_j(x)}_{\neq 0 \hbox{ only for }j=k} + \sum_{j\in\If} c_j\underbrace{\basphi_{\nu(j)}(\xno{k})}_{=0,\ k\not\in\If} = U_k $$
Some examples will further clarify the notation. With a regular left-to-right numbering of nodes in a mesh with P1 elements, and Dirichlet conditions at \( x=0 \), we use finite element basis functions associated with the nodes \( 1, 2, \ldots, N_n \), implying that \( \nu(j)=j+1 \), \( j=0,\ldots,N \), where \( N=N_n-1 \). For the particular mesh below the expansion becomes
$$ u(x) = U_0\basphi_0(x) + c_0\basphi_1(x) + c_1\basphi_2(x) + \cdots + c_4\basphi_5(x)\tp $$
Here is a mesh with an irregular cell and node numbering:
Say we in this latter mesh have Dirichlet conditions on the left-most and right-most node, with numbers 3 and 1, respectively. Then we can number the unknowns at the interior nodes from left to right, giving \( \nu(0)=0 \), \( \nu(1)=4 \), \( \nu(2)=5 \), \( \nu(3)=2 \). This gives
$$ B(x) = U_3\basphi_3(x) + U_1\basphi_1(x),$$ and
$$ u(x) = B(x) + \sum_{j=0}^3 c_j\basphi_{\nu(j)} = U_3\basphi_3 + U_1\basphi_1 + c_0\basphi_0 + c_1\basphi_4 + c_2\basphi_5 + c_3\basphi_2\tp$$
Switching to the more standard case of left-to-right numbering and boundary conditions \( u(0)=C \), \( u(L)=D \), we have \( N=N_n-2 \) and
$$ \begin{align*} u(x) &= C\basphi_0 + D\basphi_{N_n} + \sum_{j\in\If} c_j\basphi_{j+1}\\ &= C\basphi_0 + D\basphi_{N_n} + c_0\basphi_1 + c_1\basphi_2 +\cdots + c_N\basphi_{N_n-1}\tp \end{align*} $$
The idea of constructing \( B \) described here generalizes almost trivially to 2D and 3D problems: \( B=\sum_{j\in\Ifb}U_j\basphi_j \), where \( \Ifb \) is the index set containing the numbers of all the nodes on the boundaries where Dirichlet values are prescribed.
Let us see how the model problem \( -u''=2 \), \( u(0)=C \), \( u(L)=D \), is affected by a \( B(x) \) to incorporate boundary values. Inserting the expression
$$ u(x) = B(x) + \sum_{j\in\If}c_j\baspsi_j(x)$$ in \( -(u'',\baspsi_i)=(f,\baspsi_i) \) and integrating by parts results in a linear system with
$$ A_{i,j} = \int_0^L \baspsi_i'(x)\baspsi_j'(x) \dx,\quad b_i = \int_0^L (f(x) - B'(x))\baspsi_i(x) \dx\tp $$ We choose \( \baspsi_i=\basphi_{i+1} \), \( i=0,\ldots,N=N_n-2 \) if the node numbering is from left to right. (Later we also need the assumption that the cells too are numbered from left to right.) The boundary function becomes
$$ B(x) = C\basphi_0(x) + D\basphi_{N_n}(x)\tp$$ The expansion for \( u(x) \) is
$$ u(x) = B(x) + \sum_{j\in\If} c_j\basphi_{j+1}(x)\tp $$ We can write the matrix and right-hand side entries as
$$ A_{i-1,j-1} = \int_0^L \basphi_i'(x)\basphi_j'(x) \dx,\quad b_{i-1} = \int_0^L (f(x) - C\basphi_{0}'(x) - D\basphi_{N_n}'(x))\basphi_i(x) \dx, $$ for \( i,j = 1,\ldots,N+1=N_n-1 \). Note that we have here used \( B'=C\basphi_0' + D\basphi_{N_n}' \).
Most of the terms in the linear system have already been computed so we concentrate on the new contribution from the boundary function. The integral \( C\int_0^L \basphi_{0}'(x))\basphi_i(x) \dx \) can only get a nonzero contribution from the first cell, \( \Omega^{(0)}=[\xno{0},\xno{1}] \) since \( \basphi_{0}'(x)=0 \) on all other cells. Moreover, \( \basphi_{0}'(x)\basphi_i(x) \dx \neq 0 \) only for \( i=0 \) and \( i=1 \) (but \( i=0 \) is excluded), since \( \basphi_{i}=0 \) on the first cell if \( i>1 \). With a similar reasoning we realize that \( D\int_0^L \basphi_{N_n}'(x))\basphi_i(x) \dx \) can only get a nonzero contribution from the last cell. From the explanations of the calculations in the section Calculating the linear system we then find that
$$ \int_0^L \basphi_{0}'(x)\basphi_{1}(x) \dx = \frac{1}{h}\cdot\frac{1}{h} = -{\half},\quad \int_0^L \basphi_{N_n}'(x)\basphi_{N_n-1}(x) \dx = \frac{1}{h}\cdot\frac{1}{h} = \half\tp$$ The extra boundary term because of \( B(x) \) boils down to adding \( C/2 \) to \( b_{0} \) and \( -D/2 \) to \( b_{N} \).
As an equivalent alternative, we now turn to cellwise computations. The element matrices and vectors are calculated as the section Cellwise computations, so we concentrate on the impact of the new term involving \( B(x) \). We observe that \( C\basphi_0'=0 \) on all cells except \( e=0 \), and \( D\basphi_{N_n}'=0 \) on all cells except \( e=N_e \). In this case there is only one unknown in these cells since \( u(0) \) and \( u(L) \) are prescribed, so the element vector has only one entry. The entry for the last cell, \( e=N_e \), becomes
$$ \tilde b_0^{(e)} = \int_{-1}^1 \left( f - D\frac{2}{h}\frac{d\refphi_1}{dX}\right) \refphi_0\frac{h}{2} \dX = (\frac{h}{2}(2 - D\frac{2}{h}\half) \int_{-1}^1 \refphi_0 \dX = h - D/2\tp $$ Similar computations on the first cell yield
$$ \tilde b_0^{(0)} = \int_{-1}^1 \left(f - C\frac{2}{h} \frac{d\refphi_0}{dX}\right) \refphi_1\frac{h}{2} \dX = (\frac{h}{2}(2 + C\frac{2}{h}\half) \int_{-1}^1 \refphi_1 \dX = h + C/2\tp $$ When assembling these contributions, we see that \( b_0 \) gets right-hand side of the linear system gets an extra term \( C/2 \) and \( b_{N} \) gets \( -D/2 \), as in the computations in the physical domain.
From an implementational point of view, there is a convenient alternative to adding the \( B(x) \) function and using only the basis functions associated with nodes where \( u \) is truly unknown. Instead of seeking
$$ \begin{equation} u(x) = \sum_{j\in\Ifb} U_j\basphi_j(x) + \sum_{j\in\If}c_j\basphi_{\nu(j)}(x), \tag{79} \end{equation} $$ we use the sum over all degrees of freedom, including the known boundary values:
$$ \begin{equation} u(x) = \sum_{j\in\If}c_j\basphi_j(x)\tp \tag{80} \end{equation} $$ Note that the collections of unknowns \( \sequencei{c} \) in (79) and (80) are different: in (79) \( N \) counts the number of nodes where \( u \) is not known, while in (79) \( N \) counts all the nodes (\( N=N_n \)).
The idea is to compute the entries in the linear system as if no Dirichlet values are prescribed. Afterwards, we modify the linear system to ensure that the known \( c_j \) values are incorporated.
A potential problem arises for the boundary term \( [u'v]_0^L \) from the integration by parts: imagining no Dirichlet conditions means that we no longer require \( v=0 \) at Dirichlet points, and the boundary term is then nonzero at these points. However, when we modify the linear system, we will erase whatever the contribution from \( [u'v]_0^L \) should be at the Dirichlet points in the right-hand side of the linear system. We can therefore safely forget \( [u'v]_0^L \) at any point where a Dirichlet condition applies.
Let us redo the computations in the example in the section General construction of a boundary function. We solve \( -u''=2 \) with \( u(0)=0 \) and \( u(L)=D \). The expressions for \( A_{i,j} \) and \( b_i \) are the same, but the numbering is different as the numbering of unknowns and nodes now coincide:
$$ A_{i,j} = \int_0^L \basphi_i'(x)\basphi_j'(x) \dx,\quad b_{i} = \int_0^L f(x)\basphi_i(x) \dx, $$ for \( i,j = 0,\ldots,N=N_n \). The integrals involving basis functions corresponding to interior mesh nodes, \( i,j=1,\ldots,N_n-1 \), are obviously the same as before. We concentrate on the contributions from \( \basphi_0 \) and \( \basphi_{N_n} \):
$$ \begin{align*} A_{0,0} &= \int_0^L (\basphi_0')^2\dx = \int_{0}^{\xno{1}} = (\basphi_0')^2\dx \frac{1}{h},\\ A_{0,1} &= \int_0^L \basphi_0'\basphi_1'\dx = \int_{0}^{\xno{1}} \basphi_0'\basphi_1'\dx = -\frac{1}{h},\\ A_{N,N} &= \int_0^L (\basphi_0')^2\dx = \int_{\xno{N_n-1}}^{\xno{N_n}} (\basphi_0')^2\dx = \frac{1}{h},\\ A_{N,N-1} &= \int_0^L (\basphi_0')^2\dx =\int_{\xno{N_n-1}}^{\xno{N_n}} (\basphi_0')^2\dx = -\frac{1}{h}\tp \end{align*} $$
The new terms on the right-hand side are also those involving \( \basphi_0 \) and \( \basphi_{N_n} \):
$$ \begin{align*} b_0 &= \int_0^L 2\basphi_0(x) \dx = \int_0^{\xno{1}} 2\basphi_0(x)\dx = h,\\ b_N &= \int_0^L 2\basphi_{N_n}\dx = \int_{\xno{N_n-1}}^{\xno{N_n}} 2\basphi_{N_n}\dx = h\tp \end{align*} $$
The complete matrix system, involving all degrees of freedom, takes the form
$$ \begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} 1 & -1 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & -1 & 1 \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} h \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ h \end{array} \right) \tag{81} \end{equation} $$
Incorporation of Dirichlet values can now be done by replacing the first and last equation by \( c_0=0 \) and \( c_N=D \). This action changes the system to
$$ \begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} h & 0 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ -1 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h\\ D \end{array} \right) \tag{82} \end{equation} $$
Note that because we do not require \( \basphi_i(0)=0 \) and \( \basphi_i(L) \), \( i\in\If \), the boundary term \( [u'v]_0^L \) gives in principle contributions \( u'(0)\basphi_0(0) \) to \( b_0 \) and \( u'(L)\basphi_N(L) \) to \( b_N \) (\( u'\basphi_i \) vanishes for \( x=0 \) or \( x=L \) for \( i=1,\ldots,N-1 \)). Nevertheless, we erase these contributions in \( b_0 \) and \( b_N \) and insert boundary values instead. This argument shows why we can drop computing \( [u'v]_0^L \) at Dirichlet nodes when we implement the Dirichlet values by modifying the linear system.
The original matrix system (74) is symmetric, but the modifications in (82) destroy the symmetry. Our described modification will in general destroy an initial symmetry in the matrix system. This is not a particular computational disadvantage for tridiagonal systems arising in 1D problems, but may be more serious in 2D and 3D problems when the systems are large and exploiting symmetry can be important for halving the storage demands, speeding up computations, and/or making the solution algorithm more robust. Therefore, an alternative modification which preserves symmetry is frequently applied.
Let \( c_k \) be a coefficient corresponding to a known value \( u(\xno{k}) = U_k \). We want to replace equation \( k \) in the system by \( c_k=U_k \), i.e., insert zeroes in row number \( k \) in the coefficient matrix, set 1 on the diagonal, and replace \( b_k \) by \( U_k \). A symmetry-preserving modification consists in first subtracting column number \( k \) in the coefficient matrix, i.e., \( A_{i,k} \) for \( i\in\If \), times the boundary value \( U_k \), from the right-hand side: \( b_i \leftarrow b_i - A_{i,k}U_k \). Then we put zeroes in row number \( k \) and column number \( k \) in the coefficient matrix, and finally set \( b_k=U_k \). The steps in algorithmic form becomes
$$ \begin{equation} \frac{1}{h}\left( \begin{array}{ccccccccc} h & 0 & 0 &\cdots & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & 2 & -1 & \ddots & & & & & \vdots \\ 0 & -1 & 2 & -1 & \ddots & & & & \vdots \\ \vdots & \ddots & & \ddots & \ddots & 0 & & & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & \ddots & & \vdots \\ \vdots & & & 0 & -1 & 2 & -1 & \ddots & \vdots \\ \vdots & & & & \ddots & \ddots & \ddots &\ddots & 0 \\ \vdots & & & & &\ddots & \ddots &\ddots & 0 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & h \end{array} \right) \left( \begin{array}{c} c_0 \\ \vdots\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ c_{N} \end{array} \right) = \left( \begin{array}{c} 0 \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h +D/h\\ D \end{array} \right) \tag{83} \end{equation} $$
The modifications of the global linear system can alternatively be done for the element matrix and vector. (The assembled system will get the value \( n \) on the main diagonal if \( n \) elements contribute to the same unknown, but the factor \( n \) will also appear on the right-hand side and hence cancel out.)
We have, in the present computational example, the element matrix and vector (77). The modifications are needed in cells where one of the degrees of freedom is known. Here, this means the first and last cell. We compute the element matrix and vector as there are no Dirichlet conditions. The boundary term \( [u'v]_0^L \) is simply forgotten at nodes that have Dirichlet conditions because the modification of the element vector will anyway erase the contribution from the boundary term. In the first cell, local degree of freedom number 0 is known and the modification becomes
$$ \begin{equation} \tilde A^{(0)} = A = \frac{1}{h}\left(\begin{array}{rr} h & 0\\ -1 & 1 \end{array}\right),\quad \tilde b^{(0)} = \left(\begin{array}{c} 0\\ h \end{array}\right)\tp \tag{84} \end{equation} $$ In the last cell we set
$$ \begin{equation} \tilde A^{(N_e)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & -1\\ 0 & h \end{array}\right),\quad \tilde b^{(N_e)} = \left(\begin{array}{c} h\\ D \end{array}\right)\tp \tag{85} \end{equation} $$
We can also perform the symmetric modification. This operation affects only the last cell with a nonzero Dirichlet condition. The algorithm is the same as for the global linear system, resulting in
$$ \begin{equation} \tilde A^{(N-1)} = A = \frac{1}{h}\left(\begin{array}{rr} h & 0\\ 0 & 1 \end{array}\right),\quad \tilde b^{(N-1)} = \left(\begin{array}{c} h + D/h\\ D \end{array}\right)\tp \tag{86} \end{equation} $$ The reader is encouraged to assemble the element matrices and vectors and check that the result coincides with the system (83).