Finite Elements for Ordinary Differential Equations

Following are the steps involved in the FEA and the corresponding method is called the FEM (Gerald and Wheatley 2003, Bhavikatti 2005):

  • 1. Subdivide the domain (interval) [a, b] into n number of subintervals, which are called elements. The nodes of those subintervals are at x0,x1 , ..., xn and joins at x1 , x2 , ., xn_!, where x0 = a and xn = b are the end nodes. We assume xi as the nodes of the interval and number the elements from 1 to n, where element 'i' runs from x{_ ! to xi. The x{ need not be evenly spaced.
  • 2. Then we apply the Galarkin method to each element separately to interpolate (subject to the differential equation) between the end nodal values, u(x{_ x) and u(x). Here, the u's are approximations to the y(xi)'s, which are the true solution to the differential equation. The nodal values are taken as c's in our adaptation.
  • 3. Use a low-degree polynomial for u(x). Here, a first-degree polynomial has been used, although quadratics or cubics may also be used. (The development for these higher-degree polynomials is more complicated than the polynomial taken here.)
  • 4. Applying the Galarkin method to element 'i', we have a pair of equations in which the unknowns are the nodal values at the ends of element 'i', which are the c's. Performing this for each element, we obtain equations that involve all the nodal values. The nodal values are then combined to give a set of equations that we can solve for the unknown nodal values. (The process of combining the separate element equations is called assembling the system.)
  • 5. These equations are adjusted for the boundary conditions and solved to get approximations to y(x) at the nodal points. Finally, we obtain intermediate values for y(x) by linear interpolation.

Consider the following differential equation and its finite element development. It involves several steps and each step is straightforward. The differential equation that will be investigated is

subject to the boundary conditions at x = a and x = b.

Step 1

Subdivide the interval [a, b] into n elements, as discussed earlier. Focus on element 'i' that runs from xi_ j to xi. For simplicity, we call the left node L and the right node R.

Step 2

We can write u(x) for element 'i' as

The N's in Equation 5.54 are first-degree (linear) Lagrangian polynomials. These shape functions are often called hat functions due to the orientation.

The shape functions N's are nothing but the functions of x and the values vary (from unity to zero) from xL to xR. It may be noted that the c's in Equation 5.54 are independent of x.

Step 3

Applying the Galarkin method to element 'i', the residual becomes

where u(x) is substituted in place of y(x). The Galarkin method sets the integral of R weighted with each of the N's (over the length of the element) to zero.

Now expanding Equations 5.56 and 5.57, we get

Step 4

Now, transform Equation 5.58 by applying integration by parts to the first integral. In the second integral, we will take Q out from the integrand as Qint, an average value within the element. Similarly, take F outside the third integral and it will be denoted as Fint. Then Equation 5.58 becomes

In the last two terms of Equation 5.60, NL = 1 at L and is zero at R, hence Equation 5.60 can be simplified as

In a similar fashion, Equation 5.59 gives Step 5

Using Equations 5.54, 5.61 and 5.62, we get Now,


In the same manner, we also get Now,

and Step 6

Substitute the values in Step 5 and rearrange to give two linear equations of the unknown and Cr:

The pair of Equations 5.67 and 5.68 is called element equations. The same procedure is to be adapted for each element to get n such pairs.

Step 7

Furthermore, combine (assemble) all the element equations together to form a system of linear equations for the problem. Point R in element 'i' is certainly the same as point L in element 'i + 1' and the c's as c0, c1, ... , cn. We also notice that the gradient (du/dx) must be the same on either sides of the join of the elements. As such, (du/dx)x=R for element 'i' equals to (du/dx)x=L in element 'i + 1'.

Finally, the results for repeating Step 7 in the set of n + 1 equations (numbered from 0 to n) give

where the diagonal elements of [X] are where

and elements above and to the left of the diagonal in rows 1 - n are

The elements of {c} are c, i = 0 to n.

The elements of {b} are

In the preceding equations, Qintii and Fintii are values of Q and F at the midpoints of element 'i'.

Step 8

Finally, adjust the set of equations of Step 6 and incorporate the boundary conditions. There are three cases out of which two cases have been discussed here. Case 1, a Dirichlet condition is specified as - y(a) = constant [and/or y(b) = constant]. Case 2, a Neumann condition is specified as - dy/dx = constant at x=a and/or x=b. (If Q=0, we cannot have a Neumann condition at both ends, because the solution would be known only within an additive constant.) [We leave Case 3, mixed conditions, which is a modification of Case 2.]

Case 1 Dirichlet Condition

In this case, c is known at the end node. Suppose this is y(a)=Ka, then the equation in row 0 is redundant, and so we remove it from the set of equations of Step 6. In the next row, we move k10 * Ka to the right-hand side (subtracting this from the element computed in Step 6). If the condition is y(b)=Kb, we do the same but with the last and next to last equations.

Case 2 Neumann Condition

In this case, c is not known at the end node. Suppose the condition is dy/dx=Ka at x=a, we retain the equation in row 0 and substitute the given value of dy /dx into the right-hand side. If the condition is dy/dx=Kb at x=b, we do the same with the last equation.

Step 9

Solve the set of equations for the unknowns c's after adjusting, in Step 8, for the boundary conditions. The c's are approximations to y(x) at the nodes. If the intermediate values of y are needed between the nodes, we obtain them by linear interpolation.


Bhat, R. and Chakraverty, S. 2003. Numerical Analysis in Engineering. Alpha Science International Ltd., U.K.

Bhavikatti, S. S. 2005. Finite Element Analysis. New Age International, New Delhi, India. Chakraverty, S. 2008. Vibration of Plates. CRC Press, Boca Rato, FL.

Gerald, C. F. and Wheatley, P. O. 2003. Applied Numerical Analysis, 7th edn. Pearson, Noida, India. Hoffmann, K. A. and Chiang, S. T. 2000. Computational Fluid Dynamics, 4th edn. Engineering Education System, Wichita, USA.

Kreyszig, E. 2010. Advanced Engineering Mathematics. Wiley, USA.

< Prev   CONTENTS   Source   Next >