Finite Element Experiments in MATLAB

This guest article was submitted by John Coady (bio below). Would you like to submit an article? If so, please see the submission guidelines. If your article is on scientific computing, plasma modeling, or basic plasma / rarefied gas research, we are interested!

In my previous post I talked about a MATLAB implementation of the Finite Element Method and gave a few examples of it solving to Poisson and Laplace equations in 2D. The examples that I provided all used piecewise linear polynomials in the Finite Element algorithm. This time I will show how the accuracy can be improved by using higher order polynomials and conduct a number of experiments to show the improvement in the results.

We will be using Finite Element Method to solve Laplace and Poisson equations of the following form.
$$ -\nabla \cdot (k(x,y) \nabla u(x,y)) = f(x,y)$$

The Finite Element method converts this continuous equation to a discrete equation of the form
$$ KU = F$$

where lower case values like u(x,y) are continuous and upper case values like U are discrete. In the discrete Poisson equation, K is the stiffness matrix of size NxN, F is the load vector of size Nx1 and U is an Nx1 vector where N is the number of nodes in the mesh. With the Finite Element method the computer will assemble the stiffness matrix K and the load vector F as well as solve the equation KU=F to determine the best approximation to U.

The Finite Element method does this conversion from continuous to discrete by choosing a finite number of test and trial functions that are simple polynomials and approximating the exact solution by a combination of those polynomial functions. We will be using linear, quadratic and cubic polynomials as test and trial functions in the Finite Element algorithm to compute an approximate solution and compare the results to the exact solution.

Discretizing Poisson

One way the Finite Element method discretizes the continuous Poisson equation is as follows. Start with the continuous Poisson equation.
$$ -\nabla \cdot (k \nabla u) = f$$

Multiply both sides by a test function \(v\) and then integrate to give the weak form of Poisson’s equation.
$$ \iint (-\nabla \cdot (k \nabla u)) v\,{d}A = \iint fv\,dA$$

Next use integration by parts to move one of the derivatives off of u and puts it on the test function v.
$$ \iint k \nabla u\nabla v\,{d}A – ku’v]_{0}^1 = \iint fv\,dxdy$$

If the boundary conditions for u’ or v are zero then it is reduced to
$$ \iint ((k \nabla u))\nabla v\,{d}A = \iint fv\,dA$$

otherwise if \(u’ \neq 0\) then the term has to stay and it can be moved over to the right hand side.
$$ \iint ((k \nabla u))\nabla v\,{d}A = \iint fv\,dA + ku’v ]_0^1$$

Using the continuous Galerkin method to discritize \(u\) we get
$$ U = U_{1}\phi_{1} + U_{2}\phi_{2} + \ldots + U_{N}\phi_{N}$$

and the derivative
$$ U’ = U_{1}\phi_{1}’ + U_{2}\phi_{2}’ + \ldots + U_{N}\phi_{N}’$$

where \(\phi_{i}\) are trial functions. We next make the test function \(v\) equal to the trial function \( \phi\) and end up with.
$$ \iint ((k (U_{1}\phi_{1}’ + U_{2}\phi_{2}’ + \ldots + U_{N}\phi_{N}’))\phi_{i}’\,{d}A = \iint f\phi_{i}\,dA $$

In matrix form we get the discrete equation
$$ KU = F$$

where
$$ K = \iint k \nabla\phi \nabla\phi \,{d}A \\[0.3em]F = \iint f\phi\,dA$$

These integrals are evaluated in the finite element method using N-point gauss quadrature. The MATLAB code from the last post used linear polynomials for test and trial functions and showed how the stiffness matrix K and the load vector F are assembled and the linear system KU=F is solved to get give an estimate of U.

As we saw in the previous post, the finite element takes as input a mesh which can be generated by a tool like distmesh. The mesh contains triangular elements which are the pieces on which the polynomials are defined, a different polynomial in each triangle. The Finite Element algorithm will calculate a contribution to the large K and F matrices for each small triangle element matrix by integrating over all triangles in the mesh, one at a time.

P1 Linear Elements

As already mentioned, the polynomials are defined on each triangle in the mesh, which form a set of basis functions for the triangular element. Each basis function has a value of 1 at its node and a value of 0 at the other nodes giving it a pyramid shape over the triangle.
$$ \phi_{i} (x_{j},y_{j}) = \left\{ \begin{array}{l l} 1 & \quad \text{i = j}\\ 0 & \quad \text{i} \neq \text{j}\\ \end{array} \right.$$

In the case of linear finite elements, the basis functions are linear polynomials of degree 1 with the designation P1.
$$ \phi_{i} (x,y) = a_{i} + b_{i}x + c_{i}y$$

Each \(\phi_{i}\) has 3 unknown coefficients \(a_{i} b_{i} c_{i} \) so with three node points on the triangle element we can generate 3 equations for each \( \phi_{i} \).
$$ a_{1} + b_{1}x_{1} + c_{1}y_{1} = 1 \\ a_{1} + b_{1}x_{2} + c_{1}y_{2} = 0 \\ a_{1} + b_{1}x_{3} + c_{1}y_{3} = 0$$

or in matrix form for the first basis function,
$$ \begin{bmatrix} 1 & x_{1} & y_{1} \\[0.3em] 1 & x_{2} & y_{2} \\[0.3em] 1 & x_{3} & y_{3} \end{bmatrix} \begin{bmatrix} a_{1} \\[0.3em] b_{1} \\[0.3em] c_{1} \end{bmatrix} = \begin{bmatrix} 1 \\[0.3em] 0 \\[0.3em] 0 \end{bmatrix} $$

For all three basis functions in the triangle we get the equation of the form
$$ \begin{bmatrix} 1 & x_{1} & y_{1} \\[0.3em] 1 & x_{2} & y_{2} \\[0.3em] 1 & x_{3} & y_{3} \end{bmatrix} \begin{bmatrix} a_{1} & a_{2} & a_{3} \\[0.3em] b_{1} & b_{2} & b_{3} \\[0.3em] c_{1} & c_{3} & c_{3} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\[0.3em] 0 & 1 & 0 \\[0.3em] 0 & 0 & 1 \end{bmatrix} $$

This has the form
$$ M M^{-1} = I$$

so we can determine the coefficients of the basis functions just by inverting the matrix M,
$$ M = \begin{bmatrix} 1 & x_{1} & y_{1} \\[0.3em] 1 & x_{2} & y_{2} \\[0.3em] 1 & x_{3} & y_{3} \end{bmatrix} \\[0.3em] M^{-1} = \begin{bmatrix} a_{1} & a_{2} & a_{3} \\[0.3em] b_{1} & b_{2} & b_{3} \\[0.3em] c_{1} & c_{2} & c_{3} \end{bmatrix} $$

P2 Quadratic Elements

$$ \phi_{i} (x_{j},y_{j}) = \left\{ \begin{array}{l l} 1 & \quad \text{i = j}\\ 0 & \quad \text{i} \neq \text{j}\\ \end{array} \right.$$

In the case of quadratic finite elements, the basis functions are quadratic polynomials of degree 2 with the designation P2.
$$ \phi_{i} (x,y) = a_{i} + b_{i}x + c_{i}y + d_{i}x^2 + e_{i}xy + f_{i}y^2$$

There are now 6 unknown coefficients, so we need to add 3 node points to each triangle element to generate 6 equations.
$$ \begin{bmatrix} 1 & x_{1} & y_{1} & x_{1}^2 & x_{1}y_{1} & y_{1}^2 \\[0.3em] 1 & x_{2} & y_{2} & x_{2}^2 & x_{2}y_{2} & y_{2}^2 \\[0.3em] 1 & x_{3} & y_{3} & x_{3}^2 & x_{3}y_{3} & y_{3}^2 \\[0.3em] 1 & x_{4} & y_{4} & x_{4}^2 & x_{4}y_{4} & y_{4}^2 \\[0.3em] 1 & x_{5} & y_{5} & x_{5}^2 & x_{5}y_{5} & y_{5}^2 \\[0.3em] 1 & x_{6} & y_{6} & x_{6}^2 & x_{6}y_{6} & y_{6}^2 \end{bmatrix} \begin{bmatrix} a_{1} \\[0.3em] b_{1} \\[0.3em] c_{1} \\[0.3em] d_{1} \\[0.3em] e_{1} \\[0.3em] f_{1} \end{bmatrix} = \begin{bmatrix} 1 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \\[0.3em] 0 \end{bmatrix} $$

We also end up with a total of 6 basis functions, with each basis having a value 1 at its assigned node and a value of zero at all the other mesh nodes in the element.
$$ \begin{bmatrix} 1 & x_{1} & y_{1} & x_{1}^2 & x_{1}y_{1} & y_{1}^2 \\[0.3em] 1 & x_{2} & y_{2} & x_{2}^2 & x_{2}y_{2} & y_{2}^2 \\[0.3em] 1 & x_{3} & y_{3} & x_{3}^2 & x_{3}y_{3} & y_{3}^2 \\[0.3em] 1 & x_{4} & y_{4} & x_{4}^2 & x_{4}y_{4} & y_{4}^2 \\[0.3em] 1 & x_{5} & y_{5} & x_{5}^2 & x_{5}y_{5} & y_{5}^2 \\[0.3em] 1 & x_{6} & y_{6} & x_{6}^2 & x_{6}y_{6} & y_{6}^2 \end{bmatrix} \begin{bmatrix} a_{1} & a_{2} & a_{3} & a_{4} & a_{5} & a_{6} \\[0.3em] b_{1} & b_{2} & b_{3} & b_{4} & b_{5} & b_{6} \\[0.3em] c_{1} & c_{2} & c_{3} & c_{4} & c_{5} & c_{6} \\[0.3em] d_{1} & d_{2} & d_{3} & d_{4} & d_{5} & d_{6} \\[0.3em] e_{1} & e_{2} & e_{3} & e_{4} & e_{5} & e_{6} \\[0.3em] f_{1} & f_{2} & f_{3} & f_{4} & f_{5} & f_{6} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\[0.3em] 0 & 1 & 0 & 0 & 0 & 0 \\[0.3em] 0 & 0 & 1 & 0 & 0 & 0 \\[0.3em] 0 & 0 & 0 & 1 & 0 & 0 \\[0.3em] 0 & 0 & 0 & 0 & 1 & 0 \\[0.3em] 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} $$

Again this has the form
$$ M M^{-1} = I$$
so we can easily determine the coefficients of all 6 basis functions in the quadratic case by inverting the matrix M. Each column of the inverted matrix is a set of coefficients for a basis function.

P3 Cubic Elements

$$ \phi_{i} (x_{j},y_{j}) = \left\{ \begin{array}{l l} 1 & \quad \text{i = j}\\ 0 & \quad \text{i} \neq \text{j}\\ \end{array} \right.$$

In the case of cubic finite elements, the basis functions are cubic polynomials of degree 3 with the designation P3
$$ \phi_{i} (x,y) = a_{i} + b_{i}x + c_{i}y + d_{i}x^2 + e_{i}xy + f_{i}y^2 + g_{i}x^3 + h_{i}x^2y + i_{i}xy^2 + j_{i}y^3 $$

We now have 10 unknowns so we need a total of 10 node points on the triangle element to determine the coefficients of the basis functions in the same manner as the P2 case.

Experiment 1

In this experiment we will solve Poisson’s equation over a square region using Finite Elements with linear, quadratic and cubic elements and compare the results. We will solve
$$ -\nabla \cdot (k(x,y) \nabla u(x,y)) = f(x,y)$$

over a 2D square shape going from lower left coordinates (0,0) to upper right coordinates (1,1) with dirichlet boundary conditions of zero. The known solution for u(x,y) with zero boundary conditions on the square is
$$ u(x,y) = xy(1-x)(1-y)$$

We will also set the value of k(x,y) in the partial differential equation to
$$ k(x,y) = 1 + xy^2$$

taking the derivative of u with respect to x and y
$$ \dfrac{\partial u}{\partial x} = y-y^2-2xy+2xy^2 \\[0.3em] \dfrac{\partial u}{\partial y} = x-2xy-x^2+2x^2y$$

Plugging these equations into equation Poisson’s equation above and you get the following for f(x,y)
$$ f(x,y) = y^4-y^3+4y^3x-4y^4x+2y-2y^2-2x^2y+6x^2y^2+2x^3y-6x^3y^2+2x-2x^2$$

Since $$ u(x,y) $ is the known solution I can plot it below and compare it with the Finite Element approximation to \(u(x,y) \).

The Finite Element Method is designed to minimize the Energy Norm
$$ \text{Energy Norm} = (\iint k*(\nabla U^* – \nabla u^*)^2\,{d}A)^{1/2}$$

which may not result in the lowest displacement norm L2.
$$ \text{L2 Norm} = (\iint (U^* – u^*)^2\,{d}A)^{1/2}$$

Where \( U^*\) is the best solution from Finite Elements and \(u^*\) is the exact solution.

The P1, P2 and P3 Finite Element Codes in Matlab for this experiment can be found here:

• Experiment 1 Linear femcode_P1
• Experiment 1 Quadratic femcode_P2
• Experiment 1 Cubic femcode_P3

These codes used the Matlab code from the previous post as a starting point which I modified to incorporate the equations of this experiment and to implement quadratic and cubic Finite Element algorithms. I also added code to calculate the Energy Norm and L2 Norm. In the quadratic case I added code to generate 3 additional mesh points per triangular element and also used a three point gauss quadrature rule to evaluate the integral over the triangle.

For the Cubic finite element code I added 7 additional nodes per triangle element and made use of a 6 point gauss quadrature rule to evaluate the integral in calculating K and F. You can find out more about using guass quadrature to evaluate integrals over a triangle in section 6.3 of this document.

The results of experiment 1 are shown in the following table:

Refining the mesh resolution produces the following:

As can be seen from the above table the Cubic elements produced the most accurate approximation to the exact solution both in terms of the Energy Norm and the L2 Norm. Refining the mesh produced more accurate results for all polynomials. Notice that the quadratic results from the original mesh are still more accurate than the linear results on the refined mesh. The accuracy of the coarser cubic mesh results were about the same level of accuracy as than the refined quadratic mesh results.

Experiment 2

In this experiment we will solve Poisson’s equation over the unit circle and compare the results.
$$ -\nabla \cdot (k(x,y) \nabla u(x,y)) = f(x,y)$$

The exact solution for u(x,y) with zero boundary conditions on the unit circle is
$$ u(x,y) = 1 – x^2 – y^2$$

We will also set the value of k(x,y) in the partial differential equation to
$$ k(x,y) = 1$$

Taking the derivative of u with respect to x and y we get
$$ \dfrac{\partial u}{\partial x} = -2x \\[0.3em] \dfrac{\partial u}{\partial y} = -2y$$

Plugging these equations into Poisson’s equation above I get the following for f(x,y)
$$ f(x,y) = 4$$

The P1, P2 and P3 Finite Element Codes in matlab for this experiment can be found here:

• Experiment 2 Linear femcode_P1
• Experiment 2 Quadratic femcode_P2
• Experiment 2 Cubic femcode_P3

The results of experiment 2 are shown in the following images and tables.

Refining the mesh resolution produced the following:

For this experiment the Cubic elements produced the lowest Energy Norm error value of 2% followed by the quadratic at 2.89% and Linear at 9.73% for the original mesh. For the refined mesh the Energy Norm was below 1 percent for Quadratic and Cubic polynomials The Linear element resulted in lowest L2 Norm followed by the Cubic L2 norm and then the Quadratic L2 norm. This result shows that the Finite Element algorithm is designed to minimize the Energy Norm and not the displacement norm or mean square error.

In this experiment, the accuracy did not improve as much as it did in Experiment 1 when using a higher order polynomial. This is due to the fact that we are using a triangular shaped mesh to approximate the curved boundary. When you have a polygonal domain such as a rectangle or square, this can be triangulated exactly with a triangular mesh leading to highly accurate solutions with higher order polynomials. In the case of a curved domain, we can only approximate its shape with a polygonal mesh and this introduces some error in the solution. The cubic polynomial solution on the curved boundary gave only a small improvement in accuracy compared to the quadratic polynomial solution. In Experiment 1 the improvement in accuracy was much greater when going from quadratic to cubic polynomial.

Experiment 3

In this experiment we will solve Laplace’s equation over the unit circle with nonzero boundary conditions and compare the results.
$$ -\nabla \cdot (k(x,y) \nabla u(x,y)) = f(x,y)$$

The known solution is
$$ u(x,y) = 3yx^2-y^3$$

We will also set the value of k(x,y) in the partial differential equation to
$$ k(x,y) = 1$$

taking the derivative of u with respect to x and y
$$ \dfrac{\partial u}{\partial x} = 6yx \\[0.3em] \dfrac{\partial u}{\partial y} = 3x^2 – 3y^2$$

Plugging these equations into the differential equation I get the following for f(x,y)
$$f(x,y) = 0$$

which is what you expect for a Laplace equation. The boundary condition around the unit circle at radius 1.0 is set to \(sin(3\theta)\).

The P1, P2 and P3 Finite Element Codes in matlab for this experiment can be found here:

• Experiment 3 Linear femcode_P1
• Experiment 3 Quadratic femcode_P2
• Experiment 3 Cubic femcode_P3

The results of experiment 3 are shown in the following images and tables.

Refining the mesh produced the following:

In this experiment the Cubic elements produced the most accurate results both in terms of Energy Norm and L2 Norm, followed by quadratic elements and then linear elements. With the refined mesh, the Cubic polynomial produced the most accurate result in terms of Energy Norm but the Linear polynomial had the lowest L2 norm. Due to the error introduced by approximating the curved boundary by a polygonal mesh, the higher order polynomials did not improve the accuracy of the solution by any significant amount. Reference 5 in the references provided below explains that the accuracy of Finite Elements on curved boundaries can be improved by using “Isoparametric quadratic triangles” or by using a nonuniform mesh with smaller elements near the curved boundary.

Experiment 4

In this experiment we will solve Laplace’s equation over a square with nonzero boundary conditions along 2 sides and compare the results.
$$ -\nabla \cdot (k(x,y) \nabla u(x,y)) = f(x,y)$$

The known solution is
$$ u(x,y) = \sin(xy)$$

We will also set the value of k(x,y) in the partial differential equation to
$$ k(x,y) = 1$$

taking the derivative of u with respect to x and y
$$ \dfrac{\partial u}{\partial x} = y\cos(xy) \\[0.3em]\dfrac{\partial u}{\partial y} = x\cos(xy)$$

Plugging these equations into the differential equation I get the following for f(x,y)
$$ f(x,y) = y^2\sin(xy) + x^2\sin(xy)$$

The boundary condition around the square shape from (0,0) to (2,2) are as follows:
$$ u(0,y) = 0 \\[0.3em] u(x,0) = 0\\[0.3em] u(2,y) = \sin(2y)\\[0.3em] u(x,2) = \sin(2x)\\[0.3em]$$

The P1, P2 and P3 Finite Element Codes in matlab for this experiment can be found here:

• Experiment 4 Linear femcode_P1
• Experiment 4 Quadratic femcode_P2
• Experiment 4 Cubic femcode_P3

The results of experiment 4 are shown in the following images and tables.

Refining the mesh produced the following:

This Experiment produced similar results to Experiment 1. For both meshes, the Cubic elements produced the most accurate results both in terms of Energy Norm and L2 Norm, followed by quadratic elements and then linear elements. Note that the cubic polynomial from the original mesh produced more accurate results than both the linear and quadratic polynomials did in the refined mesh.

Experiment 5

In this experiment we will solve Poisson’s equation over a square with zero boundary conditions and a known solution containing an exponential function.
$$ -\nabla \cdot (k(x,y) \nabla u(x,y)) = f(x,y)$$

The known solution is
$$ u(x,y) = -16*\exp(-4*((x-1)^2 + y^2)) – 16*\exp(-4*((x+1)^2 + y^2));$$

We will also set the value of k(x,y) in the partial differential equation to
$$ k(x,y) = 1$$

taking the derivative of u with respect to x and y
$$ \dfrac{\partial u}{\partial x} = 64*(2*x-2)*\exp(-4*((x-1)^2 + y^2)) + 64*(2*x+2)*\exp(-4*((x+1)^2 + y^2)) \\[0.3em] \dfrac{\partial u}{\partial y} = 64*2*y*\exp(-4*((x-1)^2 + y^2)) + 64*2*y*\exp(-4*((x+1)^2 + y^2))$$

Plugging these equations into the differential equation I get the following for f(x,y)
$$ f(x,y) = -128*\exp(-4*((x-1)^2 + y^2)) – 128*(x-1)*(-8*(x-1))*\exp(-4*((x-1)^2 + y^2)) – 128*\exp(-4*((x-1)^2 + y^2)) – 128*y*(-4*2*y)*\exp(-4*((x-1)^2 + y^2)) -128*\exp(-4*((x+1)^2 + y^2)) – 128*(x+1)*(-8*(x+1))*\exp(-4*((x+1)^2 + y^2)) – 128*\exp(-4*((x+1)^2 + y^2)) – 128*y*(-4*2*y)*\exp(-4*((x+1)^2 + y^2))$$

This solution is plotted in 3D below.

The P1, P2 and P3 Finite Element Codes in matlab for this experiment can be found here:

• Experiment 5 Linear femcode_P1
• Experiment 5 Quadratic femcode_P2
• Experiment 5 Cubic femcode_P3

The results of experiment 5 are shown in the following images and tables.

Refining the mesh produced the following:

For this experiment, on the course mesh the Energy norm error of the linear polynomial was 40% but that of the cubic polynomial was below 1%. I added an extra table refining the mesh further. For the linear case, reducing the mesh spacing from 0.4 to 0.2 to 0.1 resulted in halving the Energy Norm for each step. For the quadratic polynomial the Energy Norm was over 4 times less for each step. For the cubic polynomial the Energy Norm was over 8 times less for each step.

Summary

• Piecewise Polynomials are used in the Finite Element method to approximate solutions to partial differential equations.
• The Finite Element method minimizes the Energy Norm when it solves Poisson’s or Laplace’s equation.
• Using higher order polynomials results in more accurate solutions in the Energy Norm sense.
• A shape with a polygonal boundary that is approximated with a polygonal mesh results in highly accurate solutions when using higher order polynomials.
• A shape with a curved boundary that is approximated with a polygonal mesh introduces errors along the boundary and the improvement in accuracy is much less pronounced when using higher order polynomials.
• Isoparametric quadratic triangles as explained in Reference 5 can be used in Finite Elements to improve the accuracy of for domains with curved boundaries.

References

1. Section 3.6 The Finite Element Method in the book “Computational Science and Engineering” by Gilbert Strang
2. Section 5.4 The Finite Element Method in the book “Introduction to Applied Mathematics” by Gilbert Strang
3. DistMesh – A Simple Mesh Generator in MATLAB
4. Finite Elements in 2D video lecture.
5. Understanding and Implementing the Finite Element Method by Mark S. Gockenback