.. !split .. _fem:deq:1D:BC:nat: Boundary conditions: specified derivative ========================================= Suppose our model problem :math:`-u''(x)=f(x)` features the boundary conditions :math:`u'(0)=C` and :math:`u(L)=D`. As already indicated in the section :ref:`fem:deq:1D:varform:ex`, 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, .. math:: \int_0^L(u''(x)+f(x)){\psi}_i(x) {\, \mathrm{d}x} = 0,\quad i\in{\mathcal{I}_s}, integrating :math:`u''{\psi}_i` by parts results in .. math:: \int_0^Lu'(x)'{\psi}_i'(x) {\, \mathrm{d}x} -(u'(L){\psi}_i(L) - u'(0){\psi}_i(0)) = \int_0^L f(x){\psi}_i(x) {\, \mathrm{d}x}, \quad i\in{\mathcal{I}_s}{\thinspace .} The first boundary term, :math:`u'(L){\psi}_i(L)`, vanishes because :math:`u(L)=D`. There are two arguments for this result, explained in detail below. The second boundary term, :math:`u'(0){\psi}_i(0)`, can be used to implement the condition :math:`u'(0)=C`, provided :math:`{\psi}_i(0)\neq 0` for some :math:`i` (but with finite elements we fortunately have :math:`{\psi}_0(0)=1`). The variational form of the differential equation then becomes .. _Eq:fem:deq:1D:BC:nat:varform: .. math:: :label: fem:deq:1D:BC:nat:varform \int_0^Lu'(x){\varphi}_i'(x) {\, \mathrm{d}x} + C{\varphi}_i(0) = \int_0^L f(x){\varphi}_i(x) {\, \mathrm{d}x},\quad i\in{\mathcal{I}_s}{\thinspace .} .. _fem:deq:1D:BC:nat:uLtest: Boundary term vanishes because of the test functions ---------------------------------------------------- At points where :math:`u` is known we may require :math:`{\psi}_i` to vanish. Here, :math:`u(L)=D` and then :math:`{\psi}_i(L)=0`, :math:`i\in{\mathcal{I}_s}`. Obviously, the boundary term :math:`u'(L){\psi}_i(L)` then vanishes. The set of basis functions :math:`\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}` contains in this case all the finite element basis functions on the mesh, expect the one that is 1 at :math:`x=L`. The basis function that is left out is used in a boundary function :math:`B(x)` instead. With a left-to-right numbering, :math:`{\psi}_i = {\varphi}_i`, :math:`i=0,\ldots,N_n-1`, and :math:`B(x)=D{\varphi}_{N_n}`: .. math:: u(x) = D{\varphi}_{N_n}(x) + \sum_{j=0}^{N=N_n-1} c_j{\varphi}_j(x){\thinspace .} Inserting this expansion for :math:`u` in the variational form :eq:`fem:deq:1D:BC:nat:varform` leads to the linear system .. _Eq:fem:deq:1D:natBC: .. math:: :label: fem:deq:1D:natBC \sum_{j=0}^{N}\left( \int_0^L {\varphi}_i'(x){\varphi}_j'(x) {\, \mathrm{d}x} \right)c_j = \int_0^L\left(f(x){\varphi}_i(x) -D{\varphi}_{N_n}'(x){\varphi}_i'(x)\right) {\, \mathrm{d}x} - C{\varphi}_i(0), for :math:`i=0,\ldots,N=N_n-1`. .. _fem:deq:1D:BC:nat:uLmod: Boundary term vanishes because of linear system modifications ------------------------------------------------------------- We may, as an alternative to the approach in the previous section, use a basis :math:`\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}` which contains all the finite element functions on the mesh: :math:`{\psi}_i={\varphi}_i`, :math:`i=0,\ldots,N_n=N`. In this case, :math:`u'(L){\psi}_i(L)=u'(L){\varphi}_i(L)\neq 0` for the :math:`i` corresponding to the boundary node at :math:`x=L` (where :math:`{\varphi}_i=1`). The number of this node is :math:`i=N_n=N` if a left-to-right numbering of nodes is utilized. However, even though :math:`u'(L){\varphi}_N(L)\neq 0`, we do not need to compute this term. For :math:`i` with an equation associated with the node :math:`x_{0}=0`. According to the section :ref:`fem:deq:1D:fem:essBC:Bfunc:modsys`, this extension consists of including :math:`A_{0,0}=1/h`, :math:`A_{0,1}=-1/h`, and :math:`b_0=h`. For :math:`i>0` we have :math:`A_{i,i}=2/h`, :math:`A_{i-1,i}=A_{i,i+1}=-1/h`. Second, we need to include the extra term :math:`-C{\varphi}_i(0)` on the right-hand side. Since all :math:`{\varphi}_i(0)=0` for :math:`i=1,\ldots,N`, this term reduces to :math:`-C{\varphi}_0(0)=-C` and affects only the first equation (:math:`i=0`). We simply add :math:`-C` to :math:`b_0` such that :math:`b_0=h - C`. Third, the boundary term :math:`-\int_0^L D{\varphi}_{N_n}(x){\varphi}_i{\, \mathrm{d}x}` must be computed. Since :math:`i=0,\ldots,N=N_n-1`, this integral can only get a nonzero contribution with :math:`i=N_n-1` over the last cell. The result becomes :math:`-Dh/6`. The resulting linear system can be summarized in the form .. _Eq:fem:deq:1D:BC:nat:Aub:system: .. math:: :label: fem:deq:1D:BC:nat:Aub:system \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){\thinspace .} Next we consider the technique where we modify the linear system to incorporate Dirichlet conditions (see the section :ref:`fem:deq:1D:BC:nat:uLmod`). Now :math:`N=N_n`. The two differences from the case above is that the :math:`-\int_0^LD{\varphi}_{N_n}{\varphi}_i{\, \mathrm{d}x}` term is left out of the right-hand side and an extra last row associated with the node :math:`x_{N_n}=L` where the Dirichlet condition applies is appended to the system. This last row is anyway replaced by the condition :math:`C_N=D` or this condition can be incorporated in a symmetric fashion. Using the simplest, former approach gives .. _Eq:fem:deq:1D:BC:nat:Aub:system:mod: .. math:: :label: fem:deq:1D:BC:nat:Aub:system:mod \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 & 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} h - C \\ 2h\\ \vdots\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ 2h \\ D \end{array} \right){\thinspace .} .. Note: show equivalence between the two systems!!! Cellwise computations (2) -------------------------- Now we compute with one element at a time, working in the reference coordinate system :math:`X\in [-1,1]`. We need to see how the :math:`u'(0)=C` condition affects the element matrix and vector. The extra term :math:`-C{\varphi}_i(0)` in the variational formulation only affects the element vector in the first cell. On the reference cell, :math:`-C{\varphi}_i(0)` is transformed to :math:`-C{\tilde{\varphi}}_r(-1)`, where :math:`r` counts local degrees of freedom. We have :math:`{\tilde{\varphi}}_0(-1)=1` and :math:`{\tilde{\varphi}}_1(-1)=0` so we are left with the contribution :math:`-C{\tilde{\varphi}}_0(-1)=-C` to :math:`\tilde b^{(0)}_0`: .. _Eq:fem:deq:1D:ex1:Ab:elm:bc:nat: .. math:: :label: fem:deq:1D:ex1:Ab:elm:bc:nat \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){\thinspace .} No other element matrices or vectors are affected by the :math:`-C{\varphi}_i(0)` boundary term. There are two alternative ways of incorporating the Dirichlet condition. Following the section :ref:`fem:deq:1D:BC:nat:uLtest`, we get a :math:`1\times 1` element matrix in the last cell and an element vector with an extra term containing :math:`D`: .. _Eq:fem:deq:1D:ex1:Ab:elm:ends: .. math:: :label: fem:deq:1D:ex1:Ab:elm:ends \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), Alternatively, we include the degree of freedom at the node with :math:`u` specified. The element matrix and vector must then be modified to constrain the :math:`\tilde c_1 = c_N` value at local node :math:`r=1`: .. _Eq:fem:deq:1D:ex1:Ab:elm:bc:nat:mod: .. math:: :label: fem:deq:1D:ex1:Ab:elm:bc:nat:mod \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){\thinspace .} .. Assemble and show that it is correct