2.3. From weak to discrete form

Finding the solution to the weak form (2.11) derived in the previous section is an infinite-dimensional problem, \(u\) is unknown at every point \(x\) in the domain and we also have an infinite number of possible weight functions \(w\). In this section we will discretize the weak form such that it can be solved with a computer. With the aim to go from a continuous problem to a discrete problem, we approximate the solution by reducing the solution space to a finite number of unknowns and at the same time restricting the number of equation to a finite set.

Shape functions

This is done by assuming a certain shape for \(u(x)\) and restricting the set of functions \(w(x)\) to all functions of that same shape.

Let’s introduce a finite set of shape functions (or interpolation functions) \(N_i(x)\) with \(i\in[1,n]\), \(n\) being the number of degrees of freedom (DOFs), and say that our approximate solution \(u^h(x)\) is a linear combination of these shape functions:

(2.12)\[ u^h(x) = \sum_{i=1}^{n}N_i(x)u_i \]

where \(u_i\) is the coefficient with which \(N_i\) is multiplied. Note that \(u_i\) is not a function of \(x\), just a coefficient defined at each DOF.

In the Finite Element method, once we have the weak form at the continuous level, we define a discretization of the domain, that is a mesh with nodes and elements, see Fig. 2.2

Fig. 2.2 Triangulation of a 2D domain (left) and discretization of a 1D domain (right).

This mesh is used to define the shape functions with certain peculiarities:

• Every shape function \(N_i\) belongs to a node with coordinates \(x_i\).

• The shape function is equal to 1 at its own node and equal to zero at all other nodes.

• The shape function \(N_i\) only takes values different from zero at the set of elements attached at to the \(i\)-th node, \(S_i\).

• The shape function \(N_i\) is defined by element-wise functions (typically polynomials)

(2.13)\[\begin{split} N_i(x) = \left\{ \begin{array}{cl} p(x), & x\in\Omega_k\ ∀ k\in S_i\\ 1, & x=x_i \\ 0, & \text{otherwise} \end{array} \right. \end{split}\]

Fig. 2.3 Interpolation of a function \(u\) using piece-wise linear polynomials (left), approximated solution arround node \(i\) (top-right) and shape functions for nodes \(i-1\), \(i\) and \(i+1\) (bottom-right).

Note that the definition of the shape functions as in (2.13) has a useful property: the coefficent \(u_i\) is equal to the solution at node \(i\): \(u^h(x_i) = u_i\).

By solving for \(u^h\) instead of for \(u\) our solution space is no longer infinitely dimensional, but has become \(n\)-dimensional. If we know the \(u_i\) values for all \(i\)’s from 1 to \(n\), using the interpolation defined in (2.12), we know \(u^h\) for every \(x\) in our domain. We are making an approximation here, it may be impossible to reconstruct the true solution \(u(x)\) exactly with the simple shape functions \(N_i(x)\), but it can be shown that if we increase the number of nodes, the true \(u(x)\) will be approximated with increased accuracy.

The following widget illustrates this concept.

Note

In the widget above, we are not solving anything yet. This is just an illustration of how increasing the resolution of the mesh allows for a better and better representation of the solution. For this plot, we are setting nodal values to the exact solution and using the shape functions to interpolate. In general, we will not know the exact solution and use the finite element method to find nodal values that give a good approximation of the PDE solution. The good nodal values that are not necessarily on top of the exact solution, but the idea that increasing the number of nodes leads to a better and better approximation does still hold.

Discrete weak form

In the Finite Element method, we substitute the discretized solution \(u^h\) into the weak form equation, as derived for our 1D-rod in equation (2.11). We are solving for \(n\) values \(u_j\). For this, we need to construct \(n\) equations. Be reminded that with the “\(\forall w\)”, we had an infinite number of equations. It turns out that it is useful to use the same shape functions \(N_i(x)\) to pick \(n\) relevant equations from this infinite set. In more technical terms, we select the same vectorial space for both the approximated solution and the test functions1. We replace the \(w(x)\) with its discretized counterpart \(w^h(x)\):

\[ w^h(x) = \sum_{i=1}^{n}N_i(x)w_i \]

In alternative notation, we can collect all shape functions in a row vector or \(1\times n\) matrix:

\[ \mathbf{N} = \left[ N_1, N_2, \ldots, N_{n} \right] \]

which gives

(2.14)\[ u^h = \mathbf{Nu} \quad\text{and}\quad w^h = \mathbf{Nw} \]

Similarly, we introduce \(\mathbf{B}\) as a matrix containing the derivatives of the shape functions

\[ \mathbf{B} = \left[ \frac{\partial N_1}{\partial x}, \frac{\partial N_2}{\partial x}, \ldots, \frac{\partial N_{n}}{\partial x} \right] \]

which gives

(2.15)\[ \frac{\partial u^h}{\partial x} = \mathbf{Bu} \quad\text{and}\quad \frac{\partial w^h}{\partial x} = \mathbf{Bw} \]

Substitution of equations (2.14) and (2.15) into equation (2.11) gives

\[ \int_{0}^{L} \mathbf{Bw}EA \mathbf{Bu}\,dx = \int_0^L \mathbf{Nw}f\,dx+\mathbf{N}(L)\mathbf{w}F \]

Note that \(\mathbf{u}\) and \(\mathbf{w}\) are not a function of \(x\) and can therefore be taken out of the integrals

(2.16)\[ \mathbf{w}^T\int_{0}^{L} \mathbf{B}^TEA \mathbf{B}\,dx\mathbf{u} = \mathbf{w}^T\int_0^L \mathbf{N}^Tf\,dx+\mathbf{w}^T\mathbf{N^T}(L)F \]

where we have made use of the the equalities \(\mathbf{Bw}=\mathbf{w}^T\mathbf{B}^T\) and \(\mathbf{Nw}=\mathbf{w}^T\mathbf{N}^T\). Finally, we can divide by \(\mathbf{w}^T\) and obtain an expression of the form

\[\mathbf{Ku} = \mathbf{f}\]

with

\[ \mathbf{K} = \int_0^L \mathbf{B}^TEA\mathbf{B}\,dx \]

and

\[ \mathbf{f} = \int_0^L \mathbf{N}^Tf\,dx + \mathbf{N}^T(L)F \]

Elimination of \(\mathbf{w}^T\)

Note that the “division by \(\mathbf{w}^T\)” that was performed above is not really a valid mathematical operation. What we are actually doing in that step is described more rigorously in this box. Equation (2.16) must hold for all \(\mathbf{w}\). If we collect all terms in equation (2.16) to the left hand side, it can be written in the form

\[ \mathbf{w}^T\mathbf{v} = 0 \]

You may recognize this as an inner product between two vectors. If the inner product \(\mathbf{w}^T\mathbf{v}\) must be equal to 0 for all \(\mathbf{w}\), it follows that

\[ \mathbf{v}=\mathbf{0} \]

In making this step we go from a scalar equation to a vector equation. Specifically, we have

\[ \mathbf{v} = \mathbf{Ku}-\mathbf{f} = \mathbf{0} \]

Now we have \(n\) equations that are linear in the \(n\) degrees of freedom in the vector \(\mathbf{u}\).

Exercise

Quiz question

---

1

Selecting the same space for both solution and test functions is known as the Bubnov-Galerkin method. In some cases is useful to select different spaces for solution and trial functions, leading to the so called Petrov-Galerkin methods, which will not be covered in this chapter.