Chapter 8
Finite elements in one dimension
Context: In this chapter, some properties of finite elements are discussed using one-dimensional examples. In particular, we show how second-order PDEs can be discretized using linear elements and how boundary conditions can be incorporated into the finite elements.
8.1 Differentiability of the Basis Functions
The simplest case of a basis of finite elements in one dimension has already been introduced in the previous chapters. The basis functions are “tent” or “hat” functions, each of which is maximal (\(=1\)) on one node and then drops linearly to the two neighboring nodes. In contrast to the Fourier basis, this basis is not differentiable (in the classical sense, not even once). This actually means that we cannot naively use this basis as a series expansion for the ansatz function to solve a second-order (or even first-order) PDE in which derivatives occur. We will now show that a weak solution can be obtained in the context of a weak formulation that does not have to meet these requirements for differentiability. This means that these basis functions can be used for second-order PDEs after all.
As an example, we continue to consider the one-dimensional Poisson equation with the residual \begin {equation} R(x) = \frac {\dif ^2 \Phi }{\dif x^2} + \frac {\rho (x)}{\varepsilon }. \label {eq:strongresidual} \end {equation} In the weighted residual method, the scalar product of the residual with a test function \(v(x)\) must vanish, \((v,R)=0\). In this integral formulation, the rules of partial integration can now be used to transfer a derivative to the test function. This yields \begin {align} (v(x), R(x)) =& \int _a^b \dif x\, v(x) \left ( \frac {\dif ^2 \Phi }{\dif x^2} + \frac {\rho (x)}{\varepsilon }\right ) \label {eq:beforepartialint} \\ =& \int _a^b \dif x\, \left [ \frac {\dif }{\dif x} \left ( v(x) \frac {\dif \Phi }{\dif x} \right ) - \frac {\dif v}{\dif x} \frac {\dif \Phi }{\dif x} \right ] + \int _a^b \dif x\, v(x) \frac {\rho (x)}{\varepsilon } \\ =& \left . v(x) \frac {\dif \Phi }{\dif x} \right |_a^b - \int _a^b \dif x\, \frac {\dif v}{\dif x} \frac {\dif \Phi }{\dif x} + \int _a^b \dif x\, v(x) \frac {\rho (x)}{\varepsilon }. \label {eq:afterpartialint} \end {align}
Equation \eqref{eq:afterpartialint} now contains no second derivative. This means that both the test function \(v(x)\) and the solution \(\Phi (x)\) need only be differentiable once for all \(x\). We can therefore use linear basis functions to discretize a second-order PDE. In the following, we will denote the simulation domain as \(\Omega \). In the one-dimensional case discussed here, \(\Omega =[a,b]\).
Note: The linear basis functions are not even simply differentiable in the classical sense, because the left- and right-sided difference quotient differs at the kinks of the function. However, the so-called weak derivative \(f'(x)\) of the function \(f(x)\) exists if \begin {equation} \int _\Omega \dif x\, v(x) f'(x) = -\int _\Omega \dif x\, v'(x) f(x) \label {eq:weakderivative} \end {equation} for arbitrary (strongly differentiable) test functions \(v(x)\). The weak derivative is not unique; for example, at the kinks of the tent functions, you can assume \(0\) (or any other value) as the weak derivative. This works because this single point makes no contribution to the integral in Eq. \eqref{eq:weakderivative}. Just as the weak derivative is not unique, the weak solution of a differential equation is not unique.
In this learning material, we will not systematically distinguish between strong and weak derivatives but implicitly assume that we are operating with the weak derivative. We often perform intuitive calculations with weak derivatives, such as taking a step function as the derivative of the finite tents and a Dirac \(\delta \)-function as the derivative of the step function. This mathematically imprecise use of derivatives is common in engineering and physics and usually leads to correct results. Spaces of weakly differentiable functions are called Sobolev spaces. Sobolev spaces are always Hilbert spaces, since a scalar product is always needed for weak differentiability, see Eq. \eqref{eq:weakderivative}.
We have thus broadened our understanding of the solution of a PDE. While for the strong solution \(R(x)\equiv 0\), Eq. \eqref{eq:strongresidual}, the function \(\Phi (x)\) must be twice strongly differentiable, for the weak solution \((v,R)\equiv 0\), Eq. \eqref{eq:afterpartialint}, this function must only be once weakly differentiable. Of course, equation \eqref{eq:afterpartialint} must be fulfilled here for all singly- differentiable test functions \(v(x)\).
Note: In the collocation method, the “strong” requirement of differentiability twice is not removed. This can be seen, for example, from the fact that the Dirac \(\delta \) as a test function for the collocation method is not even differentiable in the weak sense. That is, the test functions of the collocation method are not in the set of test functions for which \((v,R)\equiv 0\) is required in the sense of the weak solution. This is the reason why the weighted residual method, and specifically the Galerkin method, has become so widespread.
8.2 Galerkin method
In the context of the weak formulation Eq. \eqref{eq:afterpartialint}, the Galerkin method is now applied. We write again \begin {equation} \Phi (x) \approx \Phi _N(x) = \sum _{n=0}^N a_n \varphi _n(x) \end {equation} as a (finite) series expansion with our finite elements, the tent functions \(\varphi _n(x)\). We first consider the case in which the test function vanishes on the boundary of the domain \(\Omega \), i.e. at \(x=a\) and \(x=b\); thus the first term in Eq. \eqref{eq:afterpartialint} vanishes. This term vanishes, for example, for periodic domains. (However, this term becomes important again when we talk about boundary conditions for non-periodic domains.)
The Galerkin condition thus becomes \begin {equation} (\varphi _k, R_N) = - \sum _n a_n \int _\Omega \dif x\, \frac {\dif \varphi _k}{\dif x} \frac {\dif \varphi _n}{\dif x} + \frac {1}{\varepsilon } \int _\Omega \dif x\, \varphi _k(x) \rho (x) = 0 \label {eq:fegaler1d} \end {equation} with unknown \(a_n\). The integrals in Eq. \eqref{eq:fegaler1d} are merely numbers; the equation is thus a system of coupled linear equations. With \begin {equation} K_{kn} = \int _\Omega \dif x\, \frac {\dif \varphi _k}{\dif x} \frac {\dif \varphi _n}{\dif x} and f_k = \frac {1}{\varepsilon } \int _\Omega \dif x\, \varphi _k(x) \rho (x) \end {equation} we thus obtain \begin {equation} \sum _n K_{kn} a_n = f_k, \label {eq:lineq} \end {equation} or in dyadic (matrix-vector) notation \begin {equation} \t {K}\cdot \v {a} = \v {f}. \label {eq:lineqdyad} \end {equation} Thus, the differential equation is transformed into an algebraic equation that can be solved numerically. The matrix \(\t {K}\) is called the system matrix or stiffness matrix. (The latter term comes from the application of finite elements in the context of structural mechanics.) The term \(\v {f}\) is often referred to as the “right hand side” (often abbreviated as “rhs”) or load vector (again from structural mechanics).
Note: In structural mechanics, which deals with the deformation of solids, \(\t {K}\) is something like a spring constant, \(\v {f}\) a force and \(\v {a}\) the displacements of the nodes, i.e. the distance to the undeformed state. This makes Eq. \eqref{eq:lineqdyad} something like the generalization of a spring law, Hooke’s law. For this reason, the symbols \(\t {K}\) and \(\v {f}\) are traditionally used. Structural mechanics is more complicated than most other numerical applications because the area on which the constitutive equations were discretized changes due to the deformation itself. This automatically leads to non-linear equations, so-called geometric non-linearities. In this case, \(\t {K}\) itself depends on \(\v {a}\).
We can now calculate the elements of the matrix \(\t {K}\) directly. For the finite element basis, the following applies \begin {equation} \frac {\dif \varphi _n(x)}{\dif x} = \left \{ \begin {array}{ll} \frac {1}{x_n - x_{n-1}} & \text {for}\; x\in [x_{n-1},x_n]\\ -\frac {1}{x_{n+1} - x_n} & \text {for}\; x\in [x_n,x_{n+1}] \\ 0 & \text {otherwise} \end {array} \right . \label {eq:finite-element-basis-derivative} \end {equation} and thus \begin {align} K_{nn} &= \frac {1}{x_n - x_{n-1}} + \frac {1}{x_{n+1} - x_n} \\ K_{n,n+1} &= -\frac {1}{x_{n+1} - x_n} \end {align}
and \(K_{kn} = 0\) for \(|n-k|>1\). The matrix \(\t {K}\) is therefore sparse, symmetrical and almost tridiagonal.
Example: For equidistant nodes with spacing \(\Delta x = x_{n+1}-x_n\) on a periodic domain, for example, we obtain \(6\) nodes (\(N=5\)): \begin {equation} \t {K} = \frac {1}{\Delta x} \begin {pmatrix} 2 & -1 & 0 & 0 & 0 & -1 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & 0 \\ 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & -1 & 2 & -1 \\ -1 & 0 & 0 & 0 & -1 & 2 \\ \end {pmatrix} \label {eq:systemmatrix1d-periodic-not-regular} \end {equation} Note that the \(-1\)s in the upper right and lower left corners (\(K_{0N}\) and \(K_{N0}\)) appear due to periodicity. The matrix is therefore not purely tridiagonal. The right-hand side \(f_k\) depends on the specific choice of the source term, i.e. the charge density \(\rho (x)\), and looks as follows: \begin {equation} \v {f} = \begin {pmatrix} (\varphi _0(x), \rho (x))/\varepsilon \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ (\varphi _5(x), \rho (x))/\varepsilon \end {pmatrix} \label {eq:rhs1d-periodic-not-regular} \end {equation}
In the course of this learning material, it has been completely neglected so far that boundary conditions are always needed to solve differential equations. Even this periodic case is not complete. In the Fourier representation, \(n=0\) (with wave vector \(q_0=0\)) represents the mean value of the Fourier series. The solution of the Poisson equation does not specify this mean value and it must therefore be given as an additional condition.
In the context of discretization using finite elements, this manifests itself in the fact that the system matrix \(\t {K}\), as written in e.g. Eq. \eqref{eq:systemmatrix1d-periodic-not-regular}, is not regular (also invertible or nonsingular). This can be seen, for example, from the fact that the determinant of \(\t {K}\) vanishes. In other words, the linear equations are not linearly independent. In the case discussed here, the rank of the matrix, i.e. the number of linearly independent rows, is exactly one less than its dimension. We can therefore remove one of these equations (or rows) and introduce the corresponding mean value condition instead. For the periodic case, this is \begin {equation} \int _\Omega \dif x\, \Phi _N(x) = \sum _n a_n \int _\Omega \dif x\, \varphi _n(x) = \sum _n a_n \Delta x = 0. \label {eq:feaverage} \end {equation} The regular system matrix then looks like this, \begin {equation} \t {K} = \frac {1}{\Delta x} \begin {pmatrix} 2 & -1 & 0 & 0 & 0 & -1 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & 0 \\ 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & -1 & 2 & -1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ \end {pmatrix} \label {eq:systemmatrix1d-periodic-regular} \end {equation} whereby \(f_N=0\) must now also apply. The right-hand side thus becomes \begin {equation} \v {f} = \begin {pmatrix} (\varphi _0(x), \rho (x))/\varepsilon \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ 0 \end {pmatrix}. \label {eq:rhs1d-periodic-regular} \end {equation} The last row corresponds to the mean value condition Eq. \eqref{eq:feaverage}, but any row could have been replaced by this condition.
8.3 Boundary Conditions
The mean value condition of the previous example is a special case of a boundary condition. In most cases, problems are treated on finite domains \(\Omega \) in which either the function value \(\Phi (x)\) or the derivative \(\dif \Phi /\dif x\) on the boundary \(\partial \Omega \) of the domain is given. (In our one-dimensional case, \(\partial \Omega =\{a,b\}\).) The first case is called a Dirichlet boundary condition, the second case is called a Neumann boundary condition. For differential equations of higher order than two, even higher derivatives could of course occur on the boundary. Combinations of Dirichlet and Neumann boundary conditions are also possible. Here we will only discuss pure Dirichlet and pure Neumann boundary conditions.
8.3.1 Dirichlet Boundary Conditions
In the one-dimensional case discussed here, the Dirichlet boundary conditions are \begin {equation} \Phi (a) \approx \Phi _N(a) = \Phi _a \quad \text {and}\quad \Phi (b) \approx \Phi _N(b) = \Phi _b. \end {equation} These conditions lead directly to \(a_0=\Phi _a\) and \(a_N=\Phi _b\). That is, the Dirichlet conditions directly fix the corresponding coefficients of the series expansion. Note that the conditions \((\varphi _0, R)=0\) and \((\varphi _N, R)=0\) have been implicitly taken from the system of equations. However, the basis functions \(\varphi _0(x)\) and \(\varphi _N(x)\) still appear in the series expansion \(\Phi _N(x)\).
Example: In our example, the system matrix with Dirichlet boundary conditions then becomes \begin {equation} \t {K} = \frac {1}{\Delta x} \begin {pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & 0 \\ 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end {pmatrix} \label {eq:systemmatrix1d-dirichlet} \end {equation} and \(f_0=\Phi _a/\Delta x\) and \(f_N=\Phi _b/\Delta x\): \begin {equation} \v {f} = \begin {pmatrix} \Phi _a/\Delta x \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ \Phi _b/\Delta x \end {pmatrix} \label {eq:rhs1d-dirichlet} \end {equation} This matrix \(\t {K}\) is regular.
The \(\Delta x\) only appears on the right side \(f_0\) and \(f_N\) because it is a prefactor in the system matrix. In an implementation, one would push \(\Delta x\) completely to the right side and it would thus completely disappear from the first and last line.
8.3.2 Neumann boundary conditions
The Neumann boundary conditions are \begin {equation} \left .\frac {\dif \Phi }{\dif x}\right |_a \approx \left .\frac {\dif \Phi _N}{\dif x}\right |_a = \Phi ^\prime _a \quad \text {and}\quad \left .\frac {\dif \Phi }{\dif x}\right |_b \approx \left .\frac {\dif \Phi _N}{\dif x}\right |_b = \Phi ^\prime _b. \end {equation} To include these conditions in our algebraic equation, we can no longer neglect the first term in Eq. \eqref{eq:afterpartialint}. This term is not relevant for the Dirichlet boundary conditions, since the basis functions \(\varphi _1(x)\) to \(\varphi _{N-1}(x)\) vanish on the boundary \(x=a,b\). The basis functions \(\varphi _0(x)\) and \(\varphi _N(x)\) do not vanish, but in the Dirichlet case they no longer appear in our set of test functions.
However, we now have to determine how to interpret the basis functions \(\varphi _0(x)\) and \(\varphi _N(x)\), which now extend beyond the boundary of the domain. A natural interpretation is to consider only the half of the “tent” that remains in the domain.
The Galerkin approach thus leads to \begin {equation} (\varphi _k, R_N) = \left . \varphi _k(x) \frac {\dif \Phi }{\dif x} \right |_a^b - \sum _n a_n \int _\Omega \dif x\, \frac {\dif \varphi _k}{\dif x} \frac {\dif \varphi _n}{\dif x} + \frac {1}{\varepsilon }(\varphi _k(x), \rho (x)) = 0. \end {equation} The additional term on the left-hand side only plays a role at \(k=0\) and \(k=N\). We obtain \begin {equation} (\varphi _0, R_N) = - \Phi ^\prime _a - \sum _n K_{0n} a_n + \frac {1}{\varepsilon }(\varphi _0(x), \rho (x)) \end {equation} with \begin {equation} K_{00} = \frac {1}{x_1-x_0} \quad \text {and}\quad K_{01} = -\frac {1}{x_1 - x_0} \end {equation} and all other \(K_{0n}=0\). We can write this again (see Eq. \eqref{eq:lineq}) with \begin {equation} f_0 = \frac {1}{\varepsilon }(\varphi _0(x), \rho (x)) - \Phi ^\prime _a \end {equation} A corresponding set of equations applies to the right boundary with \(k=N\).
Example: In our example, the system matrix for two Neumann boundary conditions then becomes \begin {equation} \t {K} = \frac {1}{\Delta x} \begin {pmatrix} 1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & 0 \\ 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & 0 & -1 & 1 \\ \end {pmatrix}. \label {eq:systemmatrix1d-neumann} \end {equation} Note that this system matrix is slightly different from the one for Dirichlet boundary conditions, Eq. This matrix is not regular and thus the problem with two Neumann boundary conditions is indeterminate. The reason for this is the same as in the periodic case: the Neumann boundary conditions do not fix the absolute value (mean) of the potential \(\Phi \). So you either need a Dirichlet boundary condition (left or right) or again the fixing of the mean value.
The right side becomes: \begin {equation} \v {f} = \begin {pmatrix} (\varphi _0(x), \rho (x))/\varepsilon - \Phi ^\prime _a \\ (\varphi _1(x), \rho (x))/\varepsilon \\ (\varphi _2(x), \rho (x))/\varepsilon \\ (\varphi _3(x), \rho (x))/\varepsilon \\ (\varphi _4(x), \rho (x))/\varepsilon \\ (\varphi _5(x), \rho (x))/\varepsilon + \Phi ^\prime _b \end {pmatrix}. \label {eq:rhs1d-neumann} \end {equation} Note that for the evaluation of \((\varphi _0, \rho )\) and \((\varphi _5, \rho )\) only half of the corresponding tent needs to be integrated.