Ordinary Differential Equation Finite Element Solver

Extensible FEA library written in C to solve ordinary differential equations given Dirichlet boundary conditions.

Problem Solved

Numerically solve the second order differential equation given the form $ \frac{d^2 y}{d x} + a \frac{d y}{d x} + by = f(x) $, where $ f(x) $ is a forcing function. Currently the library is configured to work with linear (two-node) and quadratic (three-node) 1D line elements.

Approach / Physics

The above governing equation is discretized according to the finite element method. Performing the mathematical process, the following linear equation $ [K][x] = [F] $ is solved, where $ [K] $ represents the coefficient matrix, $ [x] $ is the variable column vector, and $ [F] $ is the constant (forcing function) column matrix.

Note that the shape of the matrices and column vectors is dependant on the number of nodes for each individual element. The equation for the coefficient matrix and constant vector components is given below:

\[ K_{ij} = \int_{-1}^1 {\phi_i'(\zeta) \phi_j'(\zeta){\cdot}\frac{1}{J} + a\phi_i(\zeta) \phi_j'(\zeta) + b\phi_i(\zeta) \phi_j(\zeta){\cdot}J \ d\zeta} \] \[ F_i = \int_{-1}^1 {f(x(\zeta)){\cdot}\phi_i(\zeta){\cdot}J \ d\zeta} \]

Where $ J $ is the Jacobian of the element, and $ \phi_n $ represents the n-th shape function for the element.

Outcome / Validation

  • The solver was validated against a known analytical solution to verify that the resulting procedure is accurate.
  • For non-analytical solutions, the solver is tested using a mesh refinement convergence study.