This study aims to address the following questions:

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.

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.

In almost all approximate methods used to determine the solution of differential and integral equations, we seek a solution in the form 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 $\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 $\Omega$ and conditions on the boundary $\Gamma$ of $\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 subjected to the boundary conditions 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 $\Omega =(0,L)$ by substituting $U_N$ into Equation (1) such that 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, 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.

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 $\Omega =(0,L)$ of the problem 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].

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: 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: 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: 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):$ 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:These functionals transform the original problem into an optimization problem, where the weak form is minimized to approximate the solution.

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:

To approach this problem, we choose the approximate solution in the form 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)$. Then, multiply Equation (21) by a weight function $w\left(x\right)$ and integrate over the domain $\Omega$$=(0,L)$. Doing this, we have Integrating the term by parts, we have Substituting Equation (22) into Equation (21) then, Applying our boundary conditions, 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: where By minimizing the quadratic functional using the relation and solving, we have 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 Differentiating with respect to the $c_i$ ’s, we have Also, ${\varphi }_i=x^i(1-x)$ satisfies boundary conditions of the differential equation.

For the choice of the approximation functions in Equation (29), 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: 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 For $N=2$, we have 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:

We recall the problem with boundary conditions, 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 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 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). where $c_1^e$ and $c_2^e$ are constants. Dividing into subdivisions with end points $x_a$ and $x_b$, we have Solving both equations simultaneously, we have Substituting $c_1^e$ and $c_1^e$ in Equations (45) and (46) into Equation (42), we have Comparing Equation (47) with Equation (5), we have We can then compute Equation (38) and Equation (39) by evaluating the integrals. We have and so on.

The coefficient or stiffness matrix is given as The coefficient matrix of the two linear finite elements with $a_e=1,c_e=-1,h_e={\displaystyle \frac{1}{2}}$ is The co-efficients $f_i^e$ are evaluated as 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 where $Q_i^e$ denote force at the nodes with $Q_1^1=0,Q_2^2=0$

Therefore, the assembled equation becomes 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$.