Basic Introduction to Finite Elements Methods. This is part of project work done for MAE 207, spring 2006. UCI. Instructor: Professor S.N. Atluri Complete course work can be found on http://12000.org/my_courses/UCI_COURSES/CREDIT_COURSES/spring_2006/spring_MAE_207/

Written by Nasser Abbasi

Introduction

This is an attempt to give a basic review of Finite Elements Methods from Mathematical point of view with examples of how it can be used to numerically solve first and second order ODE's. Currently I show how to use FEM to solve first and second order ODE. I am also working on a detailed derivation and implementation using FEM to solve the 2D Poisson's equation but this work is not yet completed.

FEM is a numerical method for solving differential equations (ordinary or partial). It can also be used to solve non-linear differential equations but I have not yet studied how this is done. FEM is a more versatile numerical method than the finite difference methods for solving differential equations as it supports more easily different types of geometry and boundary conditions, in addition the solution of the differential equation found using FEM can be used at any point in the domain and not just on the grid points as the case is with finite difference methods. On the other hand FEM is more mathematically complex method, and its implementation is not as straight forward as with finite difference methods. For examples of using FDM to solve a PDE equation this is a link to a report I wrote on this subject.

We start by considering only ordinary differential equations with constant coefficients over the $x$ domain (real line).

Consider the differential equation

MATH defined over $0\leq x\leq1$ with the boundary condition MATH. In the above, only when $y\left( x\right) $ is the exact solution, call it MATH, do we have the above identity to be true.

In other words, only when MATH we can write that MATH

Such a differential equations can be represented as an operator $L$MATH If we know the exact solution, call it MATH then we write MATH

FEM is based on the weighted residual methods (WRM) where we assume that the solution of an differential equation is the sum of weighted basis functions represented by the symbol $\phi_{i}$ in here. This is in essence is similar to Fourier series, where we represent a function as a weighted sums of series made up of the basis functions which happened to be in that case the sin and cosine functions.

So the first step in solving the differential equation is to assume that the solution, called MATH can be written as

MATH

Where $q_{j}$ are unknown coefficients (the weights) to be determined. Hence the main computational part of FEM will be focused on determining these coefficients.

When we substitute this assumed solution in the original ODE such as shown in the above example, equation (1) now becomes

MATH

Where $R\left( x\right) $ is called the differential equation residual, which is a function over $x$ and in general will not be zero due to the approximate nature of our assumed solution. Our goal is to to determine the coefficients $q_{j}$ which will make $R\left( x\right) $ the minimum over the domain of the solution.

The optimal case is for $R\left( x\right) $ to be zero over the domain. One method to be able to achieve this is by forcing $R\left( x\right) $ to meet the following requirement MATH

for all possible sets of function MATH which are also defined over the same domain. The functions MATH are linearly independent from each others. If we can make $R\left( x\right) $ satisfy the above for each one of these functions, then this implies that $R\left( x\right) $ is zero. And the solution MATH will be as close as possible to the exact solution. We will find out in FEM that the more elements we use, the closer to the exact solution we get. This property of convergence when it comes to FEM is important, but not analyzed here.

Each one of these functions MATH is called a test function (or a weight function), hence the name of this method.

In the galerkin method of FEM, the test functions are chosen from the same class of functions as the trial functions as will be illustrated below.

By making $R\left( x\right) $ satisfy the above integral equation (called the weak form of the original differential equation) for $N$ number of test functions, where $N$ is the number of the unknown coefficients $q_{i}$, then we obtain a set of $N$ algebraic equations, which we can solve for $q_{i}$.

The above is the basic outline of all methods based on the weighted residual methods. The choice of the trial basis functions, and the choice of the test functions, determine the method used. Different numerical schemes use different types of trial and test functions.

In the above, the assumed solution MATH is made up of a series of trial functions (the basis). This solution is assumed to be valid over the whole domain. This is called a global trial function. In methods such as Finite Elements and Finite volume, the domain itself is descritized, and the assumed solution is made up of a series of solutions, each of which is defined over each element resulting from the discretization process.

In addition, in FEM, the unknown coefficients, called $q_{i}$ above, have a physical meaning, they are taken as the solution values at each node. The trial functions themselves are generated by using polynomial interpolation between the nodal values. The polynomial can be linear, quadratic or cubic polynomial or higher order. Lagrangian interpolation method is normally used for this step. The order of the polynomial is determined by the number of unknowns at the nodes. For example, if our goal is to determine the displacement at the nodes, then we have 2 unknowns, one at each end of the element. Hence a linear interpolation will be sufficient in this case, since a linear polynomial $a_{0}+a_{1}x$ contain 2 unknowns, the $a_{1}$ and $a_{0}$. If in addition to the displacement, we wish to also solve for the rotation at each end of the element, hence we have a total of 4 unknowns, 2 at each end of the element, which are the displacement and the rotation. Hence the minimum interpolating polynomial needed will be a cubic polynomial MATH. In the examples below, we assume that we are only solving for displacement, hence a linear polynomial will be sufficient.

At first, we will work with global trial functions to illustrate how to use weighted residual method.

Weighted Residual method. Global trial functions.

The best way to learn how to use WRM is by working over and programming some examples.

We analyze the solution in terms of errors and the effect of changing $N$ on the result.

First example. First order ODE

Given the following ODE MATH defined over $0\leq x\leq1$ with the boundary condition MATH, we wish to solve this numerically using the WRM. This ODE has an exact solution of $y=e^{x}$.

The solution using WRM will always follow these steps.

step 1

Assume a solution that is valid over the domain $0\leq x\leq1$ to be a series solution of trial (basis) functions. We start by selecting a trial functions. The assumed solution takes the form of MATH

Where MATH is the trial function which we have to choose, and $q_{j}$ are the $N$ unknown coefficients to be determined subject to a condition which will be shown below. $\tilde{y}_{0}$ is the assumed solution which needs to be valid only at the boundary conditions. Hence in this example, since we are given that the solution must be $1$ at the initial condition $x=0$, then $\tilde{y}_{0}=1$ will satisfy this boundary condition. Hence our trial solution is

MATH

step 2

Now decide on what trial function MATH to use. For this example, we can select the trial functions to be polynomials in $x$ or trigonometric functions. Let use choose a polynomial MATH, hence our assumed solution becomes

MATH

Now we need to determine the coefficients $q_{j}$, and then our solution will be complete. This is done in the following step.

step 3

Substitute the above assumed solution into the original ODE, we obtain the residual $R\left( x\right) $MATH

$R\left( x\right) $ is the ODE residual. This is the error which will result when the assumed solution is used in place of the exact solution.

Hence from (1), we find the residual to be

MATH

Our goal now is to reduce this residual to minimum. The way we achieve this is by requiring that the residual satisfies the following integral equation

MATH

The above is a set of $N$ equations. The integration is carried over the whole domain, and MATH is a weight (test) function, which we have to also select. Depending on the numerical scheme used, the test function will assume different forms.

For the Galerkin method, we select the test function to be from the same family of functions as the trial (basis) functions. Hence in this example, let us select the test function to the following polynomial MATH

step 4

Decide on a value for $N$ and solve the set of equations generated from (3). Let us pick $N=3,$ hence $R\left( x\right) $ becomes

MATH

Substitute the above in (3) we obtain

MATH

The above generates $N$ equations to solve, there are

MATH

Now carry the integration above, we obtain the following 3 equations

MATH

Which can be written in matrix form as

MATH

The solution is MATH

Hence our assumed series solution is now complete, using the above coefficients, and from equation (1) we write

MATH

Hence

MATH

Let use compare the above solution to the exact solution $y=e^{x}$ by comparing the values of the solution over a number of points. This is done using the following small Mathematica code


p1.png

To make this more useful, we can examine how the error changes as $N$ changes. The following Mathematica code determines the solution and calculates the same table as above for $N=1\cdots5$


p2.png


p3.png




This code below plots the absolute error as $N$ changes. Notice that the number of peaks in the error plot is also $N$ which is the polynomial order (the trial solution) used to approximate the exact solution, which is to be expected.


p4.png

Second example. 4$^{th}$ order ODE.

Now we will use a more complex example and repeat the above steps. We now want to numerically solve the followingMATH

defined over $0\leq x\leq1\ $with the boundary conditions MATH. This problem is taken from Professor S.N.Atluri text book 'Methods of computer modeling in engineering and the sciences' Volume 1, page 47-50. Professor Atluri used a trigonometric functions for the trial function MATH Which already satisfies the boundary conditions. For the test function, the same function as above is used hence the test (weight) function is MATH

The book above then reduces the residual equation to a symmetric form by doing integration by parts before solving it for the coefficients. In here, we will use the unsymmetrical weak form and compare the results with those shown in the above textbook. We now start again with the same steps as we did in the above example.

step 1

Select the trial solution.

MATH

$\tilde{y}_{0}=0$ as this will satisfy the boundary conditions. Hence the trial solution is

MATH

step 2

Select trial basis function MATH. As mentioned above, we select MATH , hence the trial solution is

MATH

step 3

Substitute the above assumed solution into the original ODE, we obtain the differential equation residual $R\left( x\right) $MATH

MATH

Notice the requirement above that the trial basis functions must be 4 times differentiable, which is the case here. From above we obtain

MATH

Our goal now is to reduce this residual to minimum. The way we achieve this is by requiring that the residual satisfies the following weak form integral equation

MATH

The above is a set of $N$ equations. The integration is carried over the whole domain, and MATH is a weight (test) function, which we have to also select. As mentioned above, in this problem we select the test function to be MATH

step 4

Decide on a value for $N$ and solve the set of equations generated from (3). Let us pick $N=3,$ hence $R\left( x\right) $ becomes

MATH

Hence (3) becomes

MATH

The above generates $N$ equations to solve for the coefficients $q_{i}$

MATH

Carry the integration above and simplify and solve for $q_{i}$ we obtain the numerical solution.

This below is a Mathematica code which solves this problem for different $N$ values, and compares the error as $N$ changes. The error shown is the percentage error in the solution (approximate compared to exact) for up to $N=10$. The result below agrees well with the result in Professor Atluri textbook.


p5.png

Finite element method

Example one. First order ODE, linear interpolation

Let us first summarize what we have discussed so far.

Given a differential equation defined over domain $\Omega$, we assume its solution to be of the form MATH.

The function MATH is called the $j^{th}$ basis function. MATH is called a trial function.

The function MATH is made up of functions called the shape functions $\ N_{k\text{ }}$as they are normally called in structural mechanics books.

$q_{j}\ $ are the unknown coefficients which are determined by solving $N$ set of equations generated by setting $N$ integrals of the form MATH to zero. Where $R\left( x\right) $ is the differential equation residual and MATH is the i$^{th}$ weight function where $i=1\cdots N.$ In all what follows $N$ is the taken as the number of nodes.

In FEM, we also carry the same basic process as was described above, the differences are the following:

We divide the domain itself into a number of elements. Next, the MATH function is found by assuming the solution to be an interpolation between the nodes of the element. The solution values at the nodes are the $q_{i}$ and are of course unknown except at the boundaries as given by the problem.

We start by deciding on what interpolation between the nodes to use. We will use polynomial interpolation. Then MATH will become the interpolation function.

In addition, the coefficients $q_{j}$ represent the solution at the node $j$. These are the unknowns, which we will solve for by solving the weak form integral equation as many times as there are unknowns to solve for.

By solving for the nodal values, we can then use the interpolating function again to find the solution at any point between the nodes.

This diagram illustrates the above, using the first example given above to solve a differential equation MATH with the given boundary condition of MATH and defined over $0\leq x\leq1$


interpolate.png

Using linear interpolation, then the solution $u\left( x\right) ,$ when $x$ is located in first element, is found from MATH

But the slope of the linear interpolating line over the first element is MATH, hence the above becomes

MATH

The above is the linear interpolating polynomial. We could also have used the formula of Lagrangian interpolation to arrive at the same result.

The above is the approximate solution which is valid over the first element only. Using superscript to indicate the element number, and assuming we have equal division between nodes of length say $h$, then we write

MATH

Again, the above is valid for MATH

We now do the same for the second element

MATH

The above is valid for MATH

and finally for the 3rd element

MATH

The above is valid for MATH

This is now illustrated in the following diagram


interpolate_2.png

Since our goal is to express the global approximate solution $u$ as a series sum of basis functions each multiplied by $q_{j}$, we now rewrite each of the $u^{j}$ to allow this, as follows

MATH

Again, the above is valid for MATH. Notice the use of the following notation: Since each element will have defined on it 2 shape functions, MATH and MATH, then we use a superscript to indicate the element number. Hence for element $1$, we will write MATH and $N_{2}^{1}(x)$.

We now do the same for the second element

MATH

The above is valid for MATH

and finally for the 3rd element

MATH

The above is valid for MATH

Now we can write the global trial function as follows

MATH

Hence we see that the shape function for node 1 is MATH $\,\ $ and the shape function for node 2 is MATH and the shape function for node 3 is MATH and the shape function for the last node is MATH

Hence for the first node, the shape function is MATH for the last node MATH and the shape function for any internal node is MATH

Now we can write the approximate solution as

MATH

This completes the first part, which is expressing the global approximate solution as a sum of basis functions each multiplied by an unknowns $q$ coefficients.

Diagram below illustrates the above.


interpolate_3.png

And the diagram below illustrates the numbering used, using general numbering instead of this example which only uses 3 elements.


interpolate_4.png

Next we need to determine the residual $R\left( x\right) $ by substituting the above global solution into the original differential equation.

Let the differential equation we are solving be MATHdefined over $0\leq x\leq 1$ with the initial condition MATH.

Let $N$ be the number of nodes, then we write

MATH

Now we apply the residual weak form of the differential equation, using as the test function MATH $\ \ $

Now we write the weak form of the differential equation

MATH

Hence we obtain $N\ $equations which we can solve for the $q_{i}$

Now we evaluate these $N$ equations for this example.

First we evaluate the derivatives of the shape functions. Assume in this example we divide the domain into 3 elements, hence we have

MATH

Hence in general, if we have $N$ nodes, then we can write that for the first node, MATH and the last node MATH and for any node in the middle, MATH for MATH and MATH for MATH

Looking back at the weak form integral above, we can evaluate it as follows

MATH

For the first node only, $j=1$, we obtain

MATH

But since $\phi_{1}$ is only none zero over MATH we can simplify the above to be (knowing that MATH due to the range of integration limits, and MATH) we obtain

MATH

The above simplifies to MATH

Since $x_{2}-x_{1}=h$ and $x_{1}-x_{2}=-h$ the above becomes

MATH

The above equation gives the first row in the global stiffness matrix for any first order linear ODE of the form MATH. We see that we only need to perform numerical integration on the term MATH.

Next we determine the last equation, which will be the last row of the stiffness matrix.

For the last node only $j=N$, we obtain

MATH

Since $\phi_{N}$ domain of influence is MATH, the above simplifies to

MATH

Now since MATH the above becomes

MATH

Which simplifies to

MATH

Replace $x_{N-1}-x_{N}=-h$ the above becomes

MATH

Hence we see that the last line of the stiffness matrix can be determined directly except for the term under the integral which needs to be evaluated using numerical integration.

Now we need to determine the equation that represents any internal node. This will be any row in the global stiffness matrix between the first and the last row.

For any $j$ other than 1 or $N$ we write

MATH

Where the integral was broken into 2 parts to handle the domain of influence of the shape functions.

MATH

For MATH MATH

For MATH MATH

Hence the weak form integral can be written as MATH

which simplifies to

MATH

For equal distance between elements, MATH the simplifies to

MATH

combine we obtain

MATH

The above gives the expression for any row in the stiffness matrix other than the first and the last row.

Hence now we can write the global stiffness matrix as

MATH

Hence we see that the stiffness matrix can be build quickly without the use of any numerical integration. Integration is needed to evaluate the force (or load) vector. Depending on the forcing function, this can be simple or difficult to do.

Once the load vector is calculated, we can solve for the unknowns $q_{i}$. But before we do that, we must first replace $q_{1}$ by the initial condition given in the problem, and we must remove the first row after that as follows

MATH

Now remove the first row (but remember to multiply $u_{0}$ by the first entry in the second row

MATH

And now remove the first column after moving the first entry in the first row to the RHS to become part of the load vector

MATH

Now we solve the system above for MATH

Once the $q_{i}$ are found, we have our solution, which is

MATH

The Appendix shows a Mathematica code which solve the above general first order ODE. 4 first order ODE solved for illustration and the solution is compared to the exact solution. Animation was made to help illustrate this. the RMS error was calculated. The animations can be access by clicking on the links below

First animation.

Second animation.

Third animation.

Fourth animation.

Example Two. 2nd order ODE, Boundary value problem. Linear interpolation. Symmetric weak form.

Now we show the FEM formulation for a boundary value, second order ODE with constant coefficients of the form

MATH

with the initial conditions MATH (called essential dirichlet condition), and the neumann boundary condition MATH at $L$, for $x_{0}\leq x\leq L$

As before, we start by assuming a solution of the form MATH

And we seek to solve for the unknown $q$ by solving $N$ algebraic equations

MATH

Where $R\left( x\right) $ is the ODE residual obtained by substituting the assumed solution into the original ODE. Hence the above becomes (For simplicity, I will not keep writing $j=1\cdots N$ each time as it is assumed to be the case.

MATH

Applying integration by parts on MATH. Since MATH then integrate both sides of the above we obtain

MATH

Hence MATH

Substitute the above in equation A1 we obtain

MATH

Consider the term MATH

First consider the term MATH Since MATH at $L$ then this term becomes MATH. But MATH is non zero only for the $N^{th}$ shape function evaluated at $x=L$. Since we are using linear interpolation, the $N^{th}$ shape function is MATH which have the value of $1$ at $x=x_{N}$. Hence MATH

Now consider the term MATH, since at $x=x_{0}$ all shape functions will be zero except for $\phi _{1}$ which has the value of $1$ at $x=x_{0}$. Hence this simplifies to MATH Now recall that MATH. But at $x=x_{0}$ only $\phi _{1}$ is defined and its has the slope of $\frac{-1}{h}$, hence MATH

However,$q_{1}=u_{0}$ since that is the by definition the initial condition. Hence finally we obtainMATHhence the symmetric weak form equation A2 can now be simplified more and it becomes

MATH

We will use the trial function obtain by linear interpolation, which are shown in this table again

MATH

We know construct the global stiffness matrix.

For the first equation (which corresponds to the first row in the global stiffness matrix, we obtain)

MATH

Now since MATH, and since $\phi_{1}$ domain of influence is only from $x_{1}$ to $x_{2}$ we write the above as

MATH

But MATH and MATH and over the domain from $x_{1}$ to $x_{2},$ MATH hence the above becomes

MATH

The above simplifies to

MATH

The above gives the first row in the stiffness matrix. Since we are assuming each element will have the same length, hence $x_{1}-x_{2}=-h,$ and the above becomes

MATH

Now we obtain the last equation, which will be the last row in the global stiffness matrix

MATH

Now since MATH, and since $\phi_{N}$ domain of influence is only from $x_{N-1}$ to $x_{N}$ we write the above as

MATH

But MATH and MATH and over the domain from $x_{N-1}$ to $x_{N},$ MATH hence the above becomes

MATH

The above simplifies to

MATH

The above represents the last row in the global stiffness matrix. Since we are assuming each element will have the same length, hence $x_{N-1}-x_{N}=-h,$ and the above becomes

MATH

Now we derive the expression for any row in between the first and the last rows. For a general node $j$ we write

MATH

Break the integral into halves to make it easier to write the trial functions over the domain of influence of each we obtain

MATH

Consider the first domain MATH we have

MATH

Over this range, MATH hence MATH where MATH hence the above becomes

MATH

Now consider the second domain MATH we have

MATH

Over this range, MATH hence MATH where MATH hence the above becomes

MATH

Combine A4 and A5 and simplify we obtain

MATH

Since we are assuming each element will have the same length, hence $x_{j-1}-x_{j}=-h,$ $x_{j}-x_{j+1}=-h$ and the above becomes

MATH

The above gives any row in the stiffness matrix other than the first and the last row. Now we can write the global stiffness matrix as

MATH

we must first replace $q_{1}$ by the initial condition given in the problem, and we must remove the first row after that as follows

MATH

Multiply the first element in the second row by $u_{0}$ and remove the first row we get

MATH

Now move the first element of the first row above to the RHS and remove the first column we obtain

MATH

Now we solve for the vector MATH and this completes our solution.

A Matlab implementation is below which solves any second order ODE. This code was used to generate an animation. This animation can be accessed by clicking on the link below.

Animation

Appendix

Mathematica code to solve first order ODE using FEM

The notebook can be downloaded from here

Matlab code to solve second order ODE using FEM

The driver file is here.

The main FEM matlab function is here.

References

1. Methods of computer modeling in engineering and the sciences. Volume 1. By Professor Satya N. Atluri. Tech Science Press.

2. Class lecture notes. MAE 207. Computational methods. UCI. Spring 2006. Instructor: Professor SN Atluri.

3. Computational techniques for fluid dynamics, Volume I. C.A.J.Fletcher. Springer-Verlag

$\ \ $