$$ \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\vexd}[1]{{v_{\small\mbox{e}, #1}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\Ddt}[1]{\frac{D #1}{dt}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\Std}[1]{\hbox{Std}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \newcommand{\X}{\boldsymbol{X}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\w}{\boldsymbol{w}} \newcommand{\V}{\boldsymbol{V}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\strain}{\boldsymbol{\varepsilon}} \newcommand{\stressc}{{\sigma}} \newcommand{\strainc}{{\varepsilon}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\iz}{\boldsymbol{i}_z} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\If}{\mathcal{I}_s} % for FEM \newcommand{\Ifd}{{I_d}} % for FEM \newcommand{\Ifb}{{I_b}} % for FEM \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\sequencei}[1]{\left\{ {#1}_i \right\}_{i\in\If}} \newcommand{\basphi}{\varphi} \newcommand{\baspsi}{\psi} \newcommand{\refphi}{\tilde\basphi} \newcommand{\psib}{\boldsymbol{\psi}} \newcommand{\sinL}[1]{\sin\left((#1+1)\pi\frac{x}{L}\right)} \newcommand{\xno}[1]{x_{#1}} \newcommand{\Xno}[1]{X_{(#1)}} \newcommand{\yno}[1]{y_{#1}} \newcommand{\Yno}[1]{Y_{(#1)}} \newcommand{\xdno}[1]{\boldsymbol{x}_{#1}} \newcommand{\dX}{\, \mathrm{d}X} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integerp}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} $$

Boundary conditions: specified derivative

Suppose our model problem $ -u''(x)=f(x) $ features the boundary conditions $ u'(0)=C $ and $ u(L)=D $. As already indicated in the section Examples on variational formulations, the former condition can be incorporated through the boundary term that arises from integration by parts. This details of this method will now be illustrated in the context of finite element basis functions.

The variational formulation

Starting with the Galerkin method,

$$ \begin{equation*} \int_0^L(u''(x)+f(x))\baspsi_i(x) \dx = 0,\quad i\in\If, \end{equation*} $$ integrating $ u''\baspsi_i $ by parts results in

$$ \begin{equation*} \int_0^Lu'(x)'\baspsi_i'(x) \dx -(u'(L)\baspsi_i(L) - u'(0)\baspsi_i(0)) = \int_0^L f(x)\baspsi_i(x) \dx, \quad i\in\If\tp \end{equation*} $$

The first boundary term, $ u'(L)\baspsi_i(L) $, vanishes because $ u(L)=D $. There are two arguments for this result, explained in detail below. The second boundary term, $ u'(0)\baspsi_i(0) $, can be used to implement the condition $ u'(0)=C $, provided $ \baspsi_i(0)\neq 0 $ for some $ i $ (but with finite elements we fortunately have $ \baspsi_0(0)=1 $). The variational form of the differential equation then becomes

$$ \begin{equation*} \int_0^Lu'(x)\basphi_i'(x) \dx + C\basphi_i(0) = \int_0^L f(x)\basphi_i(x) \dx,\quad i\in\If\tp \tag{87} \end{equation*} $$

Boundary term vanishes because of the test functions

At points where $ u $ is known we may require $ \baspsi_i $ to vanish. Here, $ u(L)=D $ and then $ \baspsi_i(L)=0 $, $ i\in\If $. Obviously, the boundary term $ u'(L)\baspsi_i(L) $ then vanishes.

The set of basis functions $ \sequencei{\baspsi} $ contains in this case all the finite element basis functions on the mesh, expect the one that is 1 at $ x=L $. The basis function that is left out is used in a boundary function $ B(x) $ instead. With a left-to-right numbering, $ \baspsi_i = \basphi_i $, $ i=0,\ldots,N_n-1 $, and $ B(x)=D\basphi_{N_n} $:

$$ u(x) = D\basphi_{N_n}(x) + \sum_{j=0}^{N=N_n-1} c_j\basphi_j(x)\tp$$

Inserting this expansion for $ u $ in the variational form (87) leads to the linear system

$$ \begin{equation} \sum_{j=0}^{N}\left( \int_0^L \basphi_i'(x)\basphi_j'(x) \dx \right)c_j = \int_0^L\left(f(x)\basphi_i(x) -D\basphi_{N_n}'(x)\basphi_i(x)\right) \dx - C\basphi_i(0), \tag{88} \end{equation} $$ for $ i=0,\ldots,N=N_n-1 $.

Boundary term vanishes because of linear system modifications

We may, as an alternative to the approach in the previous section, use a basis $ \sequencei{\baspsi} $ which contains all the finite element functions on the mesh: $ \baspsi_i=\basphi_i $, $ i=0,\ldots,N_n=N $. In this case, $ u'(L)\baspsi_i(L)=u'(L)\basphi_i(L)\neq 0 $ for the $ i $ corresponding to the boundary node at $ x=L $ (where $ \basphi_i=1 $). The number of this node is $ i=N_n=N $ if a left-to-right numbering of nodes is utilized.

However, even though $ u'(L)\basphi_N(L)\neq 0 $, we do not need to compute this term. For \( iModification of the linear system or the section Symmetric modification of the linear system.

The expansion for $ u $ reads

$$ \begin{equation*} u(x) = \sum_{j\in\If} c_j\basphi_j(x), \quad B(x) = D\basphi_N(x), \end{equation*} $$ with $ N=N_n $. Insertion in the variational form (87) leads to the linear system

$$ \begin{equation} \sum_{j\in\If}\left( \int_0^L \basphi_i'(x)\basphi_j'(x) \dx \right)c_j = \int_0^L\left(f(x)\basphi_i(x)\right) \dx - C\basphi_i(0),\quad i\in\If \tp \tag{89} \end{equation} $$ After having computed the system, we replace the last row by $ c_N=D $, either straightforwardly as in Section refModification of the linear system or in a symmetric fashion as in Section refSymmetric modification of the linear system. These modifications can also be performed in the element matrix and vector for the right-most cell.

Direct computation of the global linear system

We now turn to actual computations with P1 finite elements. The focus is on how the linear system and the element matrices and vectors are modified by the condition $ u'(0)=C $.

Consider first the approach where Dirichlet conditions are incorporated by a $ B(x) $ function and the known degree of freedom $ C_{N_n} $ is left out from the linear system (see the section Boundary term vanishes because of the test functions). The relevant formula for the linear system is given by (88). There are three differences compared to the extensively computed case where $ u(0)=0 $ in the sections Computation in the global physical domain and Cellwise computations. First, because we do not have a Dirichlet condition at the left boundary, we need to extend the linear system (74) with an equation associated with the node $ \xno{0}=0 $. According to the section Modification of the linear system, this extension consists of including $ A_{0,0}=1/h $, $ A_{0,1}=-1/h $, and $ b_0=h $. For $ i>0 $ we have $ A_{i,i}=2/h $, $ A_{i-1,i}=A_{i,i+1}=-1/h $. Second, we need to include the extra term $ -C\basphi_i(0) $ on the right-hand side. Since all $ \basphi_i(0)=0 $ for $ i=1,\ldots,N $, this term reduces to $ -C\basphi_0(0)=-C $ and affects only the first equation ($ i=0 $). We simply add $ -C $ to $ b_0 $ such that $ b_0=h - C $. Third, the boundary term $ -\int_0^L D\basphi_{N_n}(x)\basphi_i\dx $ must be computed. Since $ i=0,\ldots,N=N_n-1 $, this integral can only get a nonzero contribution with $ i=N_n-1 $ over the last cell. The result becomes $ -Dh/6 $. The resulting linear system can be summarized in 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 & 2 \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 - C \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ \vdots\\ 2h - Dh/6 \end{array} \right)\tp \tag{90} \end{equation} $$

Next we consider the technique where we modify the linear system to incorporate Dirichlet conditions (see the section Boundary term vanishes because of linear system modifications). Now $ N=N_n $. The two differences from the case above is that the $ -\int_0^LD\basphi_{N_n}\basphi_i\dx $ term is left out of the right-hand side and an extra last row associated with the node $ \xno{N_n}=L $ where the Dirichlet condition applies is appended to the system. This last row is anyway replaced by the condition $ C_N=D $ or this condition can be incorporated in a symmetric fashion. Using the simplest, former approach gives

$$ \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 & -1 & 2 & -1 \\ 0 &\cdots & \cdots &\cdots & \cdots & \cdots & 0 & 0 & 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 - C \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h \\ D \end{array} \right)\tp \tag{91} \end{equation} $$

Cellwise computations

Now we compute with one element at a time, working in the reference coordinate system $ X\in [-1,1] $. We need to see how the $ u'(0)=C $ condition affects the element matrix and vector. The extra term $ -C\basphi_i(0) $ in the variational formulation only affects the element vector in the first cell. On the reference cell, $ -C\basphi_i(0) $ is transformed to $ -C\refphi_r(-1) $, where $ r $ counts local degrees of freedom. We have $ \refphi_0(-1)=1 $ and $ \refphi_1(-1)=0 $ so we are left with the contribution $ -C\refphi_0(-1)=-C $ to $ \tilde b^{(0)}_0 $:

$$ \begin{equation} \tilde A^{(0)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & 1\\ -1 & 1 \end{array}\right),\quad \tilde b^{(0)} = \left(\begin{array}{c} h - C\\ h \end{array}\right)\tp \tag{92} \end{equation} $$ No other element matrices or vectors are affected by the $ -C\basphi_i(0) $ boundary term.

There are two alternative ways of incorporating the Dirichlet condition. Following the section Boundary term vanishes because of the test functions, we get a $ 1\times 1 $ element matrix in the last cell and an element vector with an extra term containing $ D $:

$$ \begin{equation} \tilde A^{(e)} =\frac{1}{h}\left(\begin{array}{r} 1 \end{array}\right),\quad \tilde b^{(e)} = h\left(\begin{array}{c} 1 - D/6 \end{array}\right), \tag{93} \end{equation} $$

Alternatively, we include the degree of freedom at the node with $ u $ specified. The element matrix and vector must then be modified to constrain the $ \tilde c_1 = c_N $ value at local node $ r=1 $:

$$ \begin{equation} \tilde A^{(N_e)} = A = \frac{1}{h}\left(\begin{array}{rr} 1 & 1\\ 0 & 1 \end{array}\right),\quad \tilde b^{(N_e)} = \left(\begin{array}{c} h\\ D \end{array}\right)\tp \tag{94} \end{equation} $$