Advanced Solutions to Boundary Value Problems: A Comparative Study of Rayleigh-Ritz and the Finite Element Method

1. Introduction

Mathematical models are essential in science and engineering for approximating and solving real-world problems. These models, which include linear, algebraic, and differential equations, help simulate and understand complex systems based on practical observations. In structural mechanics, they are crucial for predicting material stress, strain, and deformation in design and manufacturing. Numerical methods like the Finite Element Method (FEM) have become indispensable for solving such models, particularly for problems involving complex domains, varying material properties, or nonlinear behaviors that make analytical solutions challenging [1].

The FEM provides approximate solutions by discretizing a system into smaller, manageable parts called finite elements. This subdivision allows for accurate representation of complex geometries and is computationally efficient. However, while FEM is widely applied, its comparative performance as a generalization of other methods, such as the Rayleigh-Ritz method, particularly in terms of approximation accuracy and computational efficiency, is a topic of continued research.

1.1. Research Objectives

This study aims to address the following questions:

• How do the Rayleigh-Ritz and finite element methods compare in solving boundary value problems?
• What are the sources of error in each method, and how do they affect the accuracy of the solutions?
• To what extent can the use of linear elements in FEM impact the overall performance in comparison to the Rayleigh-Ritz method?

1.2. Scope and Methodology

The study will involve solving a boundary value problem using both the Rayleigh-Ritz method and the linear element approach of the finite element method. The accuracy of the approximations will be compared through graphical illustrations and numerical analysis. Sources of error, including element approximation and domain discretization, will be identified and evaluated. While FEM will be explored using linear elements, more advanced elements such as quadratic or cubic will not be considered in this analysis. Additionally, the study will focus solely on static, linear problems and will not account for dynamic or nonlinear behavior, which could be addressed in future research.

1.3. Paper Structure

The remainder of this paper is structured as follows: Section 2 reviews the methodology of the research using the finite element method and variational principles. Section 3 outlines the mathematical formulation of both methods for the boundary value problem under consideration. Section 4 presents the comparative analysis of the results, including error estimation and graphical illustrations. Section 5 discusses the findings, and Section 6 concludes the study with suggestions for further research and practical implications.

2. Method

2.1. Weighted Integrals

In almost all approximate methods used to determine the solution of differential and integral equations, we seek a solution in the form

$$u\left(x\right)\approx U_N\left(x\right)=\sum _{j=1}^N{c}_j{\varphi }_j\left(x\right)$$

where u represents the solution of a differential equation and associated boundary conditions, and $U_N$ is its approximation that is represented as a linear combination of unknown parameters $c_j$ and known functions ${\varphi }_j$ of position x in the domain $\mathsf{\Omega }$ on which the problem is posed. We shall shortly discuss the conditions on ${\varphi }_j$. The approximate solution $U_N$ is completely known only when $c_j$ are known. Thus, we must find a means to determine $c_j$ such that $U_N$ satisfies the differential equation at every point x of the domain $\mathsf{\Omega }$ and conditions on the boundary $\mathsf{\Gamma }$ of $\mathsf{\Omega }$, then $U_N\left(x\right)=u\left(x\right)$, which is the exact solution of the problem. Of course, approximate methods are not about problems for which exact solutions can be determined by some methods of mathematical analysis: the role of approximate methods is to find an approximate solution of problems that prove difficult to obtain analytically [2].

When the exact solution cannot be determined, the alternative is to find a solution $U_N$ that satisfies the governing equations in an alternative way. In the process of satisfying the governing equations approximately, we obtain(not accidentally but by planning) N algebraic relations among the N parameters $c_1,c_2,\dots ,c_N$. For example, consider the problem of solving the differential equation

$$-{\displaystyle \frac{d}{dx}}\left[a\left(x\right){\displaystyle \frac{du}{dx}}\right]=f\left(x\right)\mathrm{for}0<x<L$$

(1)

subjected to the boundary conditions

$$u\left(0\right)=u_0,{\left[a\left(x\right){\displaystyle \frac{du}{dx}}\right]}_{x=1}=Q_0$$

where $a\left(x\right),c\left(x\right)$, and $f\left(x\right)$ are known functions, $u_0$ and $Q_0$ are known parameters, and $u\left(x\right)$ is the function to be determined. We now seek an approximate solution over the entire domain $\mathsf{\Omega }=(0,L)$ by substituting $U_N$ into Equation (1) such that

$$-{\displaystyle \frac{d}{dx}}\left[a\left(x\right){\displaystyle \frac{d{U}_N}{dx}}\right]-f\left(x\right)\equiv 0$$

(2)

We shall consider how to solve such equations later. The equation requires the approximate solution $U_N$ to satisfy the differential equation in the weighted-integral sense,

$${\int }_0^{1}w\left(x\right)Rdx=0$$

(3)

where R is called the residual defined as $R\equiv -{\displaystyle \frac{d}{dx}}\left[a\left(x\right){\displaystyle \frac{d{U}_N}{dx}}\right]-f\left(x\right)$ and $w\left(x\right)$ is called a weight function.

2.2. Development of Weak Forms

There are three steps in the development of the weak form of any differential equation.

Step 1: This step is the same as in a weighted-residual method. Move all terms of the differential equation to one side (so that it reads $\dots =0$ ), multiply the entire equation with a function $w\left(x\right)$, and integrate over the domain $\mathsf{\Omega }=(0,L)$ of the problem

$$0={\int }_0^{L}w\left[-{\displaystyle \frac{d}{dx}}\left(a{\displaystyle \frac{du}{dx}}\right)-f\right]dx$$

(4)

Recall that the expression in the square brackets is not identically zero since u is replaced by its approximation, $U_N$. Mathematically, in Equation (2), the error in the differential equation (due to the approximation of the solution) is made zero in the weighted integral sense.

Step 2: While the weighted integral statement, Equation (2), allows us to obtain the necessary number $\left(N\right)$ of algebraic relations among $c_j$ for N different choices of the weight function w, it requires that the approximation functions ${\varphi }_j$ be such that $U_N$ is differentiable as many times as needed in the original differential equation and satisfies the specified boundary conditions. So, it makes sense to shift half of the derivatives from u to w so that both are differentiated equally, and we have weaker continuity requirements on ${\varphi }_j$. The resulting integral form is known as the weak form.

Step 3: The third and last step of the weak formulation is to impose the actual boundary conditions of the problem under consideration. It is here that we require the weight function w to vanish at boundary points where the essential boundary conditions are specified, i.e., w is required to satisfy the homogeneous form of the specified essential boundary conditions of the problem [1].

2.3. Linear and Bilinear Functionals

Linear and bilinear functionals are fundamental concepts in functional analysis and variational methods, including the Rayleigh-Ritz method, which is used to approximate solutions to boundary value problems (BVPs). These functionals play a key role in formulating and solving the variational form of the governing differential equations [3]. A linear functional L maps a function u from a vector space into real numbers $\mathbb{R}$, such that:

$$L(c_1{u}_1+c_2{u}_2)=c_{1}L\left(u_1\right)+c_{2}L\left(u_2\right)$$

for any scalars $c_1,c_2$ and functions $u_1,u_2.$ In the Rayleigh-Ritz method, linear functionals are typically associated with external forces, boundary conditions, or source terms in the governing equations. A bilinear functional, B(u,v), maps two functions u and v from a vector space into real numbers $\mathbb{R}$, and it is linear in each argument:

$$B(c_1{u}_1+c_2{u}_2,v)=c_{1}B(u_1,v)+c_{2}B(u_2,v)$$

for any scalars $c_1,c_2$ and functions $u_1,u_2,andv.$ In boundary value problems, bilinear functionals often represent the energy terms in the system, such as the potential energy in elasticity or other physical quantities. The Rayleigh-Ritz method is based on variational principles, where the objective is to minimize a functional (often representing the total energy of the system) to find an approximate solution to a BVP. This is done by approximating the solution as a linear combination of trial functions and applying the variational principle. In the variational formulation of a boundary value problem, the solution satisfies a weak form of the governing differential equations. This weak form is expressed using bilinear and linear functionals as shown below:

$$\prod \left(u\right)=B(u,v)-L\left(u\right)$$

where $\prod \left(u\right)$ is the second-order differential equation for the system, B(u,v) is the bilinear functional representing the system’s internal energy, and L(u) is the linear functional representing external forces. In the Rayleigh-Ritz method, the trial solution u(x) is expressed as a linear combination of known basis functions ${\varphi }_i\left(x\right):$

$$u\left(x\right)=\sum _{i=1}^n{c}_i{\varphi }_i\left(x\right)$$

where $c_i$ are unknown coefficients to be determined. Substituting this into the variational form leads to a system of algebraic equations for $c_i$, which are obtained by minimizing the functional:

$${\displaystyle \frac{\partial \prod \left(u\right)}{\partial c_i}}=0.$$

These functionals transform the original problem into an optimization problem, where the weak form is minimized to approximate the solution.

2.4. A brief look at the Finite Element Method

The finite element method (FEM) is a numerical technique for solving problems which are described by partial differential equations. The finite element method is a technique in which a given domain is represented as a collection of simple domains, called finite elements, so that it is possible to systematically construct the approximation functions needed in a variational or weighted-residual approximation of the solution of a problem over each element [4]. Thus, the finite element method differs from the traditional Ritz, Galerkin, least-squares, collocation and other weighted residual methods in the manner in which the approximation functions are constructed. But this difference is responsible for the following three basic features of the finite element method:

• Division of whole domain into sub-domains that enable a systematic derivation of the approximation functions as well as representation of complex domains.
• Derivation of approximation functions over each element.
• Assembly of elements is based on the continuity of the solution and balance of internal fluxes; the assemblage of elements results in a numerical analog of the mathematical model of the problem being analyzed [2].

2.4.1. Discretization of the Domain

In the finite element method, the domain $\mathsf{\Omega }$ of the problem is divided into a set of subintervals i.e., line elements, called finite elements. A typical element is denoted ${\mathsf{\Omega }}_e$ and it is located between points A and B with coordinates $x_a$ and $x_b$ (i.e, of length $x_b-x_a$ ). The reason for dividing a domain into a set of sub-domains is twofold. First, domains of most systems by design are a composite of geometrically materially different parts, and the solution on these sub-domains is represented by different functions that are continuous at the interfaces of these sub-domains. Therefore, it is appropriate to seek approximation of the solution over each sub-domain. Second, approximation of the solution over each element is simpler than its approximation over the entire domain. However, the number of elements into which the total domain is divided in a problem depends mainly on the geometry of the domain and on the desired accuracy of the solution [5].

2.4.2. Derivation of Element Equations

In the finite element method, we seek an approximate solution to Equation (1) over each finite element. The polynomial approximation of the solution within a typical finite element $Q^e$ is assumed to be of the form

$$u_h^e=\sum _{j=1}^n{u}_j^e{\psi }_j^e\left(x\right)$$

(5)

where $u_j^e$ are the values of the solution $u\left(x\right)$ at the nodes of the finite element $Q^e$, and are the approximation functions over the element. Next, we develop the algebraic equations among the unknown parameters, like the Ritz and Galerkin method. The main difference here is that we work with a finite element (i.e., sub-domain) as opposed to the total domain. This step results in a matrix equation of the form $\left[K^e\right]\left\{c^e\right\}=\left\{F^e\right\}$, which is called the finite element model of the original equation [1]. The derivation of finite element equations involves the following three steps:

• Construct the weighted-residual or weak form of the differential equation.
• Obtain an approximate solution over a typical finite element.
• Derive the finite element equations by substituting the approximate solution into the weighted-residual or weak form.

To obtain the weak form, we multiply the governing differential Equation (1) with a weight function w and integrate over a typical element which results into

$$0={\int }_{x_a}^{x_b}\left(a{\displaystyle \frac{du}{dx}}{\displaystyle \frac{dw}{dx}}+cwU-wf\right)dx-{\left[wa{\displaystyle \frac{du}{dx}}\right]}_{x_a}^{x_b}$$

(6)

The last step is to identify the primary and secondary variables of the weak form. This requires us to classify the boundary conditions of each differential equation into essential (or geometric) and natural (or force) boundary conditions. The classification is made uniquely by examining the boundary term appearing in the weak form (Equation (6)),

$${\left[wa{\displaystyle \frac{du}{dx}}\right]}_{x_a}^{x_b}$$

The coefficient of the weight function w in the boundary expression is called a secondary variable. The dependent unknown u in the same form as the weight function w appearing in the boundary expression is termed a primary variable. For the model at hand, the primary variable is u while the secondary variable is $a{\displaystyle \frac{du}{dx}}$. For a typical lone element, we have four boundary conditions

$$\begin{array}{c}{u}_h^e\left(x_a\right)=u_1^e\end{array}$$

(7)

$$\begin{array}{c}{\left(-a{\displaystyle \frac{du}{dx}}\right)}_{x=x_a}=Q_1^e\end{array}$$

(8)

$$\begin{array}{c}{u}_h^e\left(x_b\right)=u_2^e\end{array}$$

(9)

$$\begin{array}{c}{\left(a{\displaystyle \frac{du}{dx}}\right)}_{x=x_b}=Q_2^e\end{array}$$

(10)

where $Q_1^e$ and $Q_2^e$ mimic the compressive force and tensile force for the axial deformation of a bar respectively. If we select $u_h^e\left(x\right)$ such that it automatically satisfies the end conditions $u_h^e\left(x_a\right)=u_1^e$ and $u_h^e\left(x_b\right)=u_2^e$, then it remains that we include the remaining conditions:

$$Q_1^e={\left(-a{\displaystyle \frac{du}{dx}}\right)}_{x=x_a},Q_2^e={\left(a{\displaystyle \frac{du}{dx}}\right)}_{x=x_b}$$

(11)

in the weak form. Using Equation (11), the weak form becomes

$$0={\int }_{x_a}^{x_b}\left(a{\displaystyle \frac{du}{dx}}{\displaystyle \frac{dw}{dx}}+cwu-wf\right)dx-w\left(x_a\right)Q_1-w\left(x_b\right)Q_2$$

(12)

The finite element model based on the weak form in Equation (12) is called the weak form Galerkin finite element model [6].

2.4.3. Assemblage of Element

The final aspect of finite element analysis is to assemble all the finite elements. In deriving the element equations, we isolated a typical element (the eth element) from the mesh and formulated the variational problem (or weak form) and developed its finite element model. To obtain the finite element equations of the total problem, we must put the elements back into their original positions. In putting the elements with their nodal degrees of freedom back into their original positions, we must require that the solution $u\left(x\right)$ is uniquely defined (i.e., u is continuous) and their source terms $Q_i^e$ are balanced at the points where elements are connected to each other. Please note, if the variable u is not continuous, we do not impose its continuity; but in the problem studied the primary variable is assumed to be continuous. The assembly of elements is carried out by imposing the following two conditions:

• If the node i of element ${\mathsf{\Omega }}^e$ is connected to the node j of element ${\mathsf{\Omega }}^f$ and node k of element ${\mathsf{\Omega }}^g$, the continuity of the primary variable u requires

$$u_i^e=u_j^f=u_k^g$$

(13)

2.

For the same three elements, the balance of secondary variables at connecting nodes requires

$$Q_i^e+Q_j^f+Q_k^g=Q_I$$

(14)

where I is the global node number assigned to the nodal point that is common to the three elements and $Q_I$ is the value of externally applied source, if any (otherwise zero).

3. Main Section

In this study, we shall consider a boundary value problem as given in the equation below:

$$-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)+c\left(x\right)u=f\left(x\right),\mathrm{for}0<x<L$$

(15)

with boundary conditions,

$${\left[EA\left(x\right){\displaystyle \frac{du}{dx}}\right]}_{x=0}=-P,\mathrm{and}{\left[EA\left(x\right){\displaystyle \frac{du}{dx}}\right]}_{x=L}=P$$

(16)

where $P>0,E,A\left(x\right),c\left(x\right)$ and $f\left(x\right)$ are given data where $E=$ Young’s Modulus, $A=$ Cross Sectional Area.

3.1. Variational Method for Solving Boundary Value Problems (Ritz Method)

To approach this problem, we choose the approximate solution in the form

$$U_2=c_1{\varphi }_1+c_2{\varphi }_2+{\varphi }_0$$

(17)

with ${\varphi }_0=1,{\varphi }_1\left(x\right)=x^2-2x,{\varphi }_2\left(x\right)=x^3-3x$. Then, we construct the weak form by moving all the terms to only one side of the differential equation such that we have $(\dots =0)$.

$$0=-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)+c\left(x\right)u-f\left(x\right),\mathrm{for}0<x<L$$

(18)

Then, multiply Equation (21) by a weight function $w\left(x\right)$ and integrate over the domain $\mathsf{\Omega }$$=(0,L)$. Doing this, we have

$$\begin{array}{c}0={\int }_0^{L}w\left(\left[-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)\right]+c\left(x\right)u-f\left(x\right)\right)dx\end{array}$$

(19)

$$\begin{array}{c}0={\int }_0^L\left(w\left[-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)\right]+wc\left(x\right)u-wf\left(x\right)dx\right)\end{array}$$

(20)

$$\begin{array}{c}0={\int }_0^{L}w\left[-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)\right]dx+{\int }_0^{L}wc\left(x\right)udx-\int wf\left(x\right)dx\end{array}$$

(21)

Integrating the term

$${\int }_0^{L}w\left[-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)\right]dx$$

by parts, we have

$${\int }_0^{L}w\left[-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)\right]dx={\int }_0^{L}EA{\displaystyle \frac{du}{dx}}\left({\displaystyle \frac{dw}{dx}}\right)dx-{\left[wEA{\displaystyle \frac{du}{dx}}\right]}_0^L$$

(22)

Substituting Equation (22) into Equation (21) then,

$$0={\int }_0^L\left[EA{\displaystyle \frac{du}{dx}}{\displaystyle \frac{dw}{dx}}+wcu-wf\right]dx-{\left[wEA{\displaystyle \frac{du}{dx}}\right]}_0^L$$

(23)

Applying our boundary conditions,

$$0={\int }_0^L\left[EA{\displaystyle \frac{du}{dx}}{\displaystyle \frac{dw}{dx}}+wcu\right]dx-{\int }_0^{L}wfdx-{\left[wEA{\displaystyle \frac{du}{dx}}\right]}_0^L$$

(24)

Equation (24) is called the weak form of the equation. The word weak refers to the weakened continuity of u, which is required to be twice differentiable in the weighted integral statement Equation (21) but only once differentiable in Equation (24).

The variational problem and quadratic functional can be expressed in the form:

$$B(w,u)=I\left(w\right)$$

(25)

where

$$\begin{array}{cc}& B(w,u)={\int }_0^L\left(EA{\displaystyle \frac{dw}{dx}}{\displaystyle \frac{du}{dx}}+cwU\right)dx\hfill \\ & l\left(w\right)=-{\int }_0^{L}wfdx+(w\left(L\right)+w\left(0\right))P\hfill \end{array}$$

By minimizing the quadratic functional using the relation

$$I\left(u\right)={\displaystyle \frac{1}{2}}B(u,u)-l\left(u\right)$$

and solving, we have

$$I={\displaystyle \frac{1}{2}}{\int }_0^L\left(EA{\left({\displaystyle \frac{du}{dx}}\right)}^2+c{u}^2\right)dx-\left({\int }_0^{L}ufdx+(u\left(L\right)+u\left(0\right))\right)P$$

(26)

Setting $EA=1,c=-1,f=-x^2,L=1,P=0$. let $U\left(1\right)=0$ and substituting $u\approx U_N$ into Equation (26), we have

$$I\left(c_1,c_2,\dots ,c_N\right)={\displaystyle \frac{1}{2}}{\int }_0^1\left[{\left(\sum _{j=1}^N{c}_j{\displaystyle \frac{d{\varphi }_j}{dx}}\right)}^2-{\left(\sum _{j=1}^N{c}_j{\varphi }_j\right)}^2+\left(2x^2\right)\sum _{j=1}^N{c}_j{\varphi }_j\right]dx$$

(27)

Differentiating with respect to the $c_i$ ’s, we have

$$\begin{array}{c}{\displaystyle \frac{\partial I}{\partial c_i}}=\left[\left({\int }_0^1\left({\displaystyle \frac{d{\varphi }_i}{dx}}{\displaystyle \frac{d{\varphi }_j}{dx}}\right)-{\varphi }_i{\varphi }_j\right)dx\right]c_j+{\int }_0^1{\varphi }_i{x}^{2}dx\end{array}$$

(28)

$$\begin{array}{c}=\sum _{j=1}^N{K}_{ij}{C}_j-F_i\end{array}$$

(29)

Also, ${\varphi }_i=x^i(1-x)$ satisfies boundary conditions of the differential equation.

For the choice of the approximation functions in Equation (28), the matrix coefficients $K_{ij}=B\left({\varphi }_i,{\varphi }_j\right)$ and vector coefficients $F_i=l\left({\varphi }_i\right)$ can be computed as follows:

$$\begin{array}{cc}\hfill & K_{ij}={\int }_0^1\left(\left[i{x}^{i-1}-(i+1)x^i\right]\left[j{x}^{j-1}-(j+1)x^j\right]-\left(x^i-x^{i+1}\right)\left(x^j-x^{j-1}\right)\right)dx\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{2ij}{(i+j){(i+j)}^2-1}}-{\displaystyle \frac{2}{(i+j+1)(i+j+2)(i+j+3)}}.\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & F_i=-{\int }_0^1{x}^2\left(x^i-x^{i+1}\right)dx\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & =-{\displaystyle \frac{1}{(i+3)(i+4)}},\hfill \end{array}$$

(33)

for $i,j=1,2,\dots ,N$. We shall consider the one, two parameter approximations.

For $N=1$, we have $K_{11}={\displaystyle \frac{3}{10}},F_1=-{\displaystyle \frac{1}{20}}$ and $c_1=-{\displaystyle \frac{1}{6}}$.

The one-parameter Rayleigh-Ritz solution us given by

$$U_1=c_1{\varphi }_1=-{\displaystyle \frac{1}{6}}\left(x-x^2\right)$$

For $N=2$, we have

$$\left[\begin{array}{cc}0.3& 0.15\\ 0.15& 0.124\end{array}\right]\left[\begin{array}{c}{c}_1\hfill \\ c_2\hfill \end{array}\right]=\left[\begin{array}{c}-0.05\\ -0.03333\end{array}\right]$$

Solving the linear equation using Crammer’s rule, we obtain $c_1=-0.08197,c_2=-0.16939$ The two parameter Ritz solution is given by:

$$\begin{array}{c}{U}_2=c_1{\varphi }_1+c_2{\varphi }_2=-0.08197\left(x-x^2\right)-0.16939\left(x^2-x^3\right)\end{array}$$

(34)

$$\begin{array}{c}=-0.08197x-0.08742x^2+0.16939x^3\end{array}$$

(35)

3.2. Using the Finite Element Analysis

We recall the problem

$$-{\displaystyle \frac{d}{dx}}\left(EA\left(x\right){\displaystyle \frac{du}{dx}}\right)+c\left(x\right)u=f\left(x\right)\mathrm{for}0<x<L$$

(36)

with boundary conditions,

$${\left[EA\left(x\right){\displaystyle \frac{du}{dx}}\right]}_{x=0}=-P\mathrm{and}{\left[EA\left(x\right){\displaystyle \frac{du}{dx}}\right]}_{x=L}=P$$

where $P>0,E,A\left(x\right),c\left(x\right)$ and $f\left(x\right)$ are given data.

We divide the domain over two subdivisions i.e. we have two finite elements. Using the Ritz method, we obtain the weak form over each subdivision which is given as

$$\begin{array}{c}0={\int }_{x_a}^{x_b}\left(EA{\displaystyle \frac{du}{dx}}{\displaystyle \frac{dw}{dx}}+cwu-wf\right)dx-\left(w\left(x_b\right)+w\left(x_a\right)\right)P\end{array}$$

(37)

$$\begin{array}{c}{K}_{ij}^e=\left[{\int }_{x_a}^{x_b}EA{\displaystyle \frac{d{\psi }_i^e}{dx}}{\displaystyle \frac{d{\psi }_j^e}{dx}}+c{\psi }_i^e{\psi }_j^e\right]dx\end{array}$$

(38)

$$\begin{array}{c}{f}_i^e={\int }_{x_a}^{x_b}f{\psi }_i^{e}dx\end{array}$$

(39)

Following our previous assumption, we set $EA=1,c=-1,f=-x^2,L=1,P=0$. The coefficient matrix over a finite element is given as

$$\begin{array}{c}\left.K_{ij}^e=\left[{\int }_{x_a}^{x_b}{\displaystyle \frac{d{\psi }_i^e}{\overrightarrow{dx}}}{\displaystyle \frac{d{\psi }_j^e}{\overrightarrow{dx}}}\right)-{\psi }_i^e{\psi }_j^e\right]dx\end{array}$$

(40)

$$\begin{array}{c}{f}_i^e={\int }_{x_a}^{x_b}\left(-x^2{\psi }_i^e\right)dx\end{array}$$

(41)

Since the weak form over an element is equivalent to the differential equation and the natural boundary conditions, the approximate solution $u_h^e$ of (Equation (5)) is required to satisfy only the end conditions $u_h^e\left(x_a\right)=u_1^e$ and $u_h^e\left(x_b\right)=u_2^e$. We seek the approximate solution

in the form of algebraic polynomials. For the weak form in Equation (37), the minimum polynomial order of $u_h^e$ is linear (which is what we are going to employ to solve this problem).

$$u_h^e\left(x\right)=c_1^e+c_2^e\left(x\right)$$

(42)

where $c_1^e$ and $c_2^e$ are constants. Dividing into subdivisions with end points $x_a$ and $x_b$, we have

$$\begin{array}{cc}\hfill & u_h^e\left(x_a\right)=c_1^e+c_2^e\left(x_a\right)\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill & u_h^e\left(x_b\right)=c_1^e+c_2^e\left(x_b\right)\hfill \end{array}$$

(44)

Solving both equations simultaneously, we have

$$\begin{array}{c}{c}_1^e={\displaystyle \frac{u_1^e{x}_b-u_2^e{x}_a}{x_b-x_a}}\end{array}$$

(45)

$$\begin{array}{c}{c}_2^e={\displaystyle \frac{u_2^e-u_1^e}{x_b-x_a}}\end{array}$$

(46)

Substituting $c_1^e$ and $c_1^e$ in Equations (45) and (46) into Equation (42), we have

$$u_h^e\left(x\right)=u_1^e\left({\displaystyle \frac{x_b-x}{x_b-x_a}}\right)+u_2^e\left({\displaystyle \frac{x-x_a}{x_b-x_a}}\right)$$

(47)

Comparing Equation (47) with Equation (5), we have

$${\psi }_1^e\left(x\right)={\displaystyle \frac{x_b-x}{x_b-x_a}},{\psi }_2^e\left(x\right)={\displaystyle \frac{x-x_a}{x_b-x_a}}$$

(48)

We can then compute Equations (38) and (39) by evaluating the integrals. We have

$$\begin{array}{c}{K}_{11}^e={\displaystyle \frac{1}{h_e}}-{\displaystyle \frac{1}{3}}{h}_e\\ K_{12}^e=-{\displaystyle \frac{1}{h_e}}-{\displaystyle \frac{1}{6}}{h}_e=K_{21}^e\\ K_{22}^e=-{\displaystyle \frac{7}{3h_e}}-{\displaystyle \frac{1}{3}}{h}_e\end{array}$$

and so on.

The coefficient or stiffness matrix is given as

$$\left[K^e\right]={\displaystyle \frac{a_e}{h_e}}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]+{\displaystyle \frac{c_e{h}_e}{6}}\left[\begin{array}{cc}2\hfill & 1\hfill \\ 1\hfill & 2\hfill \end{array}\right]$$

The coefficient matrix of the two linear finite elements with $a_e=1,c_e=-1,h_e={\displaystyle \frac{1}{2}}$ is

$$\left[K^e\right]={\displaystyle \frac{1}{12}}\left[\begin{array}{cc}22& -25\\ -25& 22\end{array}\right]$$

The co-efficients $f_i^e$ are evaluated as

$$\begin{array}{cc}\hfill f_i^e& =-{\displaystyle \frac{1}{h_e}}\left[{\displaystyle \frac{x_b}{3}}\left(x_b^3-x_a^3\right)-{\displaystyle \frac{1}{4}}\left(x_b^4-x_a^4\right)\right]\hfill \\ \hfill f_2^e& =-{\displaystyle \frac{1}{h_e}}\left[{\displaystyle \frac{1}{4}}\left(x_b^4-x_a^4\right)-{\displaystyle \frac{x_a}{3}}\left(x_b^3-x_a^3\right)\right]\hfill \end{array}$$

Element $1\left(h_1={\displaystyle \frac{1}{2}},x_a=0,x_b={\displaystyle \frac{1}{2}}\right)$

$f_1^1=-0.0104167,f_2^1=-0.03125$.

Element $2\left(h_2={\displaystyle \frac{1}{2}},x_a={\displaystyle \frac{1}{2}},x_b=1\right)$

$f_1^2=-0.11458333,f_2^2=-0.177083333$.

The assembled equation is given as

$$\left[\begin{array}{ccc}{K}_{11}^1& K_{12}^1& 0\\ K_{21}^1& K_{22}^1+K_{11}^2& K_{12}^2\\ 0& K_{21}^2& K_{22}^2\end{array}\right]\left[\begin{array}{c}{u}_1^e\hfill \\ u_2^e\hfill \\ u_3^e\hfill \end{array}\right]=\left[\begin{array}{c}{f}_1^1+Q_1^1\\ \left(f_2^1+f_1^2\right)+\left(Q_2^1+Q_1^2\right)\\ f_2^2+Q_2^2\end{array}\right]$$

where $Q_i^e$ denote force at the nodes with $Q_1^1=0,Q_2^2=0$

Therefore, the assembled equation becomes

$$\left[\begin{array}{ccc}{K}_{11}^1& K_{12}^1& 0\\ K_{21}^1& K_{22}^1+K_{11}^2& K_{12}^2\\ 0& K_{21}^2& K_{22}^2\end{array}\right]\left[\begin{array}{c}{u}_1^e\hfill \\ u_2^e\hfill \\ u_3^e\hfill \end{array}\right]=\left[\begin{array}{c}{f}_1^1\\ f_2^1+f_1^2\\ f_2^2\end{array}\right]$$

$$\left[\begin{array}{ccc}1.8333& -2.0833& 0\\ -2.0833& 2.1666& -2.0833\\ 0& -2.0833& 1.8333\end{array}\right]\left[\begin{array}{c}{U}_1\hfill \\ U_2\hfill \\ U_3\hfill \end{array}\right]=\left[\begin{array}{c}-0.0104167\\ -0.145833\\ -0.1770833\end{array}\right]$$

According to the boundary condition, $U_1=0.0$ and $U_3=0.0$. Therefore, solving for $U_2$, we have the solution: $U_1=0.0,U_2=-0.06731,U_3=0.0$.

4. Results

Figure 1. The Graphical Representation of Rayleigh-Ritz Solution

Figure 1. The Graphical Representation of Rayleigh-Ritz Solution

Figure 2. Graph of the Finite Element Method in Two Subdivisions

Figure 2. Graph of the Finite Element Method in Two Subdivisions

5. Discussion

A systematic study of the steps involved in the finite element formulation of a model second order differential equation in single variable was presented. The study introduces the basic principles of the finite element method and applied to problems in one dimension. Taking a close look at Table 1, we discover that the values obtained using the finite element method is reasonably close to those obtained using the Ritz method. This shows that the finite element method is equivalent to the use of the Rayleigh-Ritz method with a piecewise polynomial approximation to the displacement. The approximation is defined through a finite set of nodal values that constitute the degrees of freedom of the problem.

Table 1. Comparison between the solutions of Rayleigh-Ritz and Finite Element Methods.

ine	Ritz Solution		FEM
ine x	$\mathrm{N}=1$	$\mathrm{N}=2$	Linear Element
ine 0.0	0.00	0.00	0.00
ine 0.1	-0.00150	-0.0089	-0.0135
ine 0.2	-0.0267	-0.0185	-0.0269
ine 0.3	-0.0350	-0.0279	-0.0404
ine 0.4	-0.0400	-0.0359	-0.0538
ine 0.5	-0.0417	-0.0417	-0.0673
ine 0.6	-0.0400	-0.0441	-0.0538
ine 0.7	-0.0350	-0.0421	-0.0404
ine 0.8	-0.0267	-0.0348	-0.0269
ine 0.9	-0.0150	-0.0211	-0.0135
ine 1.0	0.00	0.00	0.0
ine

Table 1. Comparison between the solutions of Rayleigh-Ritz and Finite Element Methods.

ine	Ritz Solution		FEM
ine x	$\mathrm{N}=1$	$\mathrm{N}=2$	Linear Element
ine 0.0	0.00	0.00	0.00
ine 0.1	-0.00150	-0.0089	-0.0135
ine 0.2	-0.0267	-0.0185	-0.0269
ine 0.3	-0.0350	-0.0279	-0.0404
ine 0.4	-0.0400	-0.0359	-0.0538
ine 0.5	-0.0417	-0.0417	-0.0673
ine 0.6	-0.0400	-0.0441	-0.0538
ine 0.7	-0.0350	-0.0421	-0.0404
ine 0.8	-0.0267	-0.0348	-0.0269
ine 0.9	-0.0150	-0.0211	-0.0135
ine 1.0	0.00	0.00	0.0
ine

6. Conclusions

In using the finite element method, the use of quadratic element is highly recommended for higher studies compared to the use of linear element for a better and more accurate solution. Also, in this research work, the domain was divided into two subdivisions to avoid too much computation. It is recommended that the domain be split into as many subdivisions as possible.