up
PDF

Examples solving an ODE using ﬁnite elements method. Mathematica and Matlab implementation and animation

Sept 28,2006 page compiled on August 19, 2015 at 10:35am

1 Introduction

This is a basic review of Finite Elements Methods from Mathematical point of view with examples of how it can be used to numerically solve ﬁrst and second order ODE’s. Currently I show how to use FEM to solve ﬁrst 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 diﬀerential equations (ordinary or partial). It can also be used to solve non-linear diﬀerential equations but I have not yet studied how this is done. FEM is a more versatile numerical method than the ﬁnite diﬀerence methods for solving diﬀerential equations as it supports more easily diﬀerent types of geometry and boundary conditions, in addition the solution of the diﬀerential 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 ﬁnite diﬀerence methods. On the other hand FEM is more mathematically complex method, and its implementation is not as straight forward as with ﬁnite diﬀerence methods.

Considering only ordinary diﬀerential equations with constant coeﬃcients over the domain (real line).

deﬁned over with the boundary condition . In the above, only when is the exact solution, call it , do we have the above identity to be true.

In other words, only when we can write that . Such a diﬀerential equations can be represented as an operator

 (1)

If we know the exact solution, call it then we write

FEM is based on the weighted residual methods (WRM) where we assume that the solution of an diﬀerential equation is the sum of weighted basis functions represented by the symbol 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 ﬁrst step in solving the diﬀerential equation is to assume that the solution, called can be written as

Where are unknown coeﬃcients (the weights) to be determined. Hence the main computational part of FEM will be focused on determining these coeﬃcients.

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

Where is called the diﬀerential equation residual, which is a function over and in general will not be zero due to the approximate nature of our assumed solution. Our goal is to to determine the coeﬃcients which will make the minimum over the domain of the solution.

The optimal case is for to be zero over the domain. One method to be able to achieve this is by forcing to meet the following requirement

for all possible sets of function which are also deﬁned over the same domain. The functions are linearly independent from each others. If we can make satisfy the above for each one of these functions, then this implies that is zero. And the solution will be as close as possible to the exact solution. We will ﬁnd 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 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 satisfy the above integral equation (called the weak form of the original diﬀerential equation) for number of test functions, where is the number of the unknown coeﬃcients , then we obtain a set of algebraic equations, which we can solve for .

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. Diﬀerent numerical schemes use diﬀerent types of trial and test functions.

In the above, the assumed solution 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 deﬁned over each element resulting from the discretization process.

In addition, in FEM, the unknown coeﬃcients, called 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 suﬃcient in this case, since a linear polynomial contain 2 unknowns, the and . 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 . In the examples below, we assume that we are only solving for displacement, hence a linear polynomial will be suﬃcient.

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

2 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 eﬀect of changing on the result.

2.1 First example. First order ODE

Given the following ODE

deﬁned over with the boundary condition , we wish to solve this numerically using the WRM. This ODE has an exact solution of .

The solution using WRM will always follow these steps.

STEP 1

Assume a solution that is valid over the domain to be a series solution of trial (basis) functions. We start by selecting a trial functions. The assumed solution takes the form of

Where is the trial function which we have to choose, and are the unknown coeﬃcients to be determined subject to a condition which will be shown below. 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 at the initial condition , then will satisfy this boundary condition. Hence our trial solution is

STEP 2

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

 (1)

Now we need to determine the coeﬃcients , 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

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 ﬁnd the residual to be

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

 (3)

The above is a set of equations. The integration is carried over the whole domain, and is a weight (test) function, which we have to also select. Depending on the numerical scheme used, the test function will assume diﬀerent 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

 (4)

STEP 4

Decide on a value for and solve the set of equations generated from (3). Let us pick hence becomes

Substitute the above in (3) we obtain

The above generates equations to solve, there are

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

Which can be written in matrix form as

The solution is

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

Hence

Let use compare the above solution to the exact solution by comparing the values of the solution over a number of points. This is done using the following small Mathematica code

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

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

2.2 Second example. 4th order ODE

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

deﬁned over with the boundary conditions . 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

Which already satisﬁes the boundary conditions. For the test function, the same function as above is used hence the test (weight) function is

The book above then reduces the residual equation to a symmetric form by doing integration by parts before solving it for the coeﬃcients. 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.

as this will satisfy the boundary conditions. Hence the trial solution is

step 2

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

step 3

Substitute the above assumed solution into the original ODE, we obtain the diﬀerential equation residual

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

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

 (3)

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

 (4)

step 4

Decide on a value for and solve the set of equations generated from (3). Let us pick hence becomes

Hence (3) becomes

The above generates equations to solve for the coeﬃcients

Carry the integration above and simplify and solve for we obtain the numerical solution.

This below is a Mathematica code which solves this problem for diﬀerent values, and compares the error as changes. The error shown is the percentage error in the solution (approximate compared to exact) for up to . The result below agrees well with the result in Professor’s Atluri textbook.

3 Finite element method

3.1 Example one. First order ODE, linear interpolation

Let us ﬁrst summarize what we have discussed so far. Given a diﬀerential equation deﬁned over domain , we assume its solution to be of the form .

The function is called the basis function. is called a trial function.

The function is made up of functions called the shape functions as they are normally called in structural mechanics books.

are the unknown coeﬃcients which are determined by solving set of equations generated by setting integrals of the form to zero. Where is the diﬀerential equation residual and is the i weight function where In all what follows is the taken as the number of nodes.

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

We divide the domain itself into a number of elements. Next, the 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 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 will become the interpolation function.

In addition, the coeﬃcients represent the solution at the node . 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 ﬁnd the solution at any point between the nodes.

This diagram illustrates the above, using the ﬁrst example given above to solve a diﬀerential equation with the given boundary condition of and deﬁned over

Using linear interpolation, then the solution when is located in ﬁrst element, is found from

But the slope of the linear interpolating line over the ﬁrst element is , hence the above becomes

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 ﬁrst element only. Using superscript to indicate the element number, and assuming we have equal division between nodes of length say , then we write

Again, the above is valid for

We now do the same for the second element

The above is valid for

and ﬁnally for the 3rd element

The above is valid for

This is now illustrated in the following diagram

Since our goal is to express the global approximate solution as a series sum of basis functions each multiplied by , we now rewrite each of the to allow this, as follows

Again, the above is valid for . Notice the use of the following notation: Since each element will have deﬁned on it 2 shape functions, and , then we use a superscript to indicate the element number. Hence for element , we will write and .

We now do the same for the second element

The above is valid for

and ﬁnally for the 3rd element

The above is valid for

The global trial function is

The shape function for node 1 is

And the shape function for node 2 is

The shape function for node 3 is

The shape function for the last node is

Hence for the ﬁrst node, the shape function is for the last node and the shape function for any internal node is

The approximate solution is

This completes the ﬁrst part, which was to express the global approximate solution as sum of basis functions, each multiplied by an unknowns coeﬃcients.

The diagram below illustrates the above.

The diagram below illustrates the numbering used.

Next the residual is found by substituting the above global solution into the original diﬀerential equation.

The diﬀerential equation to solved is

Deﬁned over with the initial condition .

is the number of nodes, therefore the residual is

The test functions are

The weak form of the diﬀerential equation is

equations are obtained from the above, which are solved for the

The derivatives of the shape functions are found ﬁrst. Assuming in this example that the domain is divided into 3 elements results in

In general, if there are nodes, then for the ﬁrst node

And for the last node

And for any node in the middle

for and

for .

Looking back at the weak form integral above, it is evaluated as follows

For the ﬁrst node only, , the following results

Since is none zero only over , and given that due to the range of integration limits, and that , then the above can simpliﬁed to give

The above simpliﬁes further to

Since and the above reduces

The above equation gives the ﬁrst row in the global stiﬀness matrix for any ﬁrst order linear ODE of the form . The above shows that numerical integration is only needed to be performed on the term

Next the last equation is found, which will be the last row of the stiﬀness matrix. For the last node only the following results

Since domain of inﬂuence is , the above simpliﬁes to

Since the above becomes

Which simpliﬁes to

Letting in the above becomes

Hence the last line of the stiﬀness matrix can be determined directly except for the term under the integral which needs to be evaluated using numerical integration.

Now the equation that represents any internal node is found. This will be any row in the global stiﬀness matrix between the ﬁrst and the last row.

For any other than 1 or the following results

Where the integral was broken into two parts to handle the domain of inﬂuence of the shape functions.

For

For

Hence the weak form integral can be written as

Which simpliﬁes to

For equal distance between elements, the above simpliﬁes to

Combining gives

The above gives the expression for any row in the stiﬀness matrix other than the ﬁrst and the last row.

Hence the global stiﬀness matrix is

The above shows that the stiﬀness matrix can be build quickly without the use of any numerical integration in this example. Integration is needed to evaluate the force (or load) vector. Depending on the forcing function, this can be simple or diﬃcult to perform.

Once the load vector is calculated, the unknowns are solved for. But before doing that, is ﬁrst replaced by the initial condition given in the problem, and the ﬁrst row is removed as follows

Now the ﬁrst row is removed (and remembering to multiply by the ﬁrst entry in the second row) gives

The ﬁrst column is now removed after moving the ﬁrst entry in the ﬁrst row to the RHS to become part of the load vector

Now the system above is solved for . Once the are found, the solution is

The Appendix shows a Mathematica code which solve the above general ﬁrst order ODE. 4 ﬁrst 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 (shows only in the HTML version of this report).

?? ?? ?? ??

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

Now the FEM formulation is given for a boundary value, second order ODE with constant coeﬃcients of the form

with the initial conditions (called essential Dirichlet condition), and the Neumann boundary condition at , for

As before, the solution is assumed to be of the form

And the unknown are found by solving algebraic equations

Where is the ODE residual obtained by substituting the assumed solution into the original ODE. Hence the above becomes (For simplicity, is not written each time as it is assumed to be the case)

Applying integration by parts on . Since

then integrating both sides of the above gives

Hence

Substituting the above in equation A1 gives

 (A2)

Consider the term

First consider the term . Since at then this term becomes . But is non zero only for the shape function evaluated at . Since linear interpolation is used, the shape function is which have the value of at . Hence

Now consider the term , since at all shape functions will be zero except for which has the value of at . Hence this simpliﬁes to

recalling that .

But at only is deﬁned and its has the slope of , hence

However, since that is by deﬁnition the initial condition. Finally one obtains

Hence the symmetric weak form equation A2 can now be simpliﬁed more giving

 (A3)

The trial function obtain by linear interpolation are used. These are shown again in this table

The global stiﬀness matrix is now constructed. For the ﬁrst equation (which corresponds to the ﬁrst row in the global stiﬀness matrix) the result is

Since , and since domain of inﬂuence is only from to then above becomes

But and and over the domain from to , hence the above becomes

The above simpliﬁes to

The above gives the ﬁrst row in the stiﬀness matrix. Since it is assumed that each element will have the same length, hence and the above becomes

For the last equation, which will be the last row in the global stiﬀness matrix

Since , and since domain of inﬂuence is only from to The above becomes

But and and over the domain from to hence the above becomes

The above simpliﬁes to

The above represents the last row in the global stiﬀness matrix. Since it is assumed that each element will have the same length, hence and the above becomes

The expression for any row in between the ﬁrst and the last rows is now found. For a general node the result is

Breaking the integral into halves to make it easier to write the trial functions over the domain of inﬂuence gives

Considering the ﬁrst domain gives

Over this range, hence where hence the above becomes

 (A4)

Now considering the second domain gives

Over this range, hence where hence the above becomes

 (A5)

Combine A4 and A5 and simplify we obtain

Since we are assuming each element will have the same length, hence and the above becomes

The above gives any row in the stiﬀness matrix other than the ﬁrst and the last row. Now we can write the global stiﬀness matrix as

we must ﬁrst replace by the initial condition given in the problem, and we must remove the ﬁrst row after that as follows

Multiply the ﬁrst element in the second row by and remove the ﬁrst row we get

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

Now we solve for the vector 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.

4 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 ﬂuid dynamics, Volume I. C.A.J.Fletcher. Springer-Verlag