# 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*x*, ...,_{0},x_{1}*x*and joins at_{n}*x*, x_{1}_{2}, .,*x*_{n}_!, where x_{0}=*a*and*x*=_{n}*b*are the end nodes. We assume*x*as the nodes of the interval and number the elements from 1 to n, where element_{i}*'i'*runs from*x*_{{}_ ! to x_{i}. 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(x_{i})'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 x_{i}_ j to x_{i}. 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 *x _{L}* to

*x*It may be noted that the c's in Equation 5.54 are independent of

_{R}.*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 *Q _{int},* an average value within the element. Similarly, take

*F*outside the third integral and it will be denoted as

*F*Then Equation 5.58 becomes

_{int}.

In the last two terms of Equation 5.60, *N _{L}* = 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,

and

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 c_{0}, *c _{1},* ... , c

_{n}. 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}=

*for element*

_{R}*'i'*equals to (du/dx)

_{x}=

*in element 'i + 1'.*

_{L}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, *Q*_{intii} and *F*_{intii} 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)*=K_{a}, 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 *k _{10}* *

*K*to the right-hand side (subtracting this from the element computed in Step 6). If the condition is y(b)=

_{a}*K*we do the same but with the last and next to last equations.

_{b},Case 2 Neumann Condition

In this case, *c* is not known at the end node. Suppose the condition is *dy/dx*=*K _{a}* 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*=

*K*at x=b, we do the same with the last equation.

_{b}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.

## Bibliography

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.