## Approximate Methods: The Rayleigh Ritz Method: Bars Under Axial Loads

#### Learning Outcomes

the steps required to find an approximate solution for a beam system using the Rayleigh Ritz method. (Step1: Assume a displacement function, apply the BC. Step 2: Write the expression for the PE of the system. Step 3: Find the minimizers of the PE of the system.)__Describe__the RR method to__Employ__an approximate solution for the displacement in a beam under axial load.__compute__between the requirement for an approximate solution and an exact solution.__Differentiate__

The Rayleigh Ritz method is a classical approximate method to find the displacement function of an object such that the it is in equilibrium with the externally applied loads. It is regarded as an ancestor of the widely used Finite Element Method (FEM). The Rayleigh Ritz method relies on the principle of minimum potential energy for conservative systems. The method involves assuming a form or a shape for the unknown displacement functions, and thus, the displacement functions would have a few unknown parameters. These assumed shape functions are termed Trial Functions. Afterwards, the potential energy function of the system is written in terms of those few parameters, and the values of those parameters that would minimize the potential energy function of the system are calculated. Mathematically speaking, we assume that the unknown displacement function is a member of a certain space of functions; for example, we can assume that where is the set of all the possible vector valued linear functions defined on the body represented by the set . With this assumption, we restrict the possible solutions to those in the form , and thus, the unknown parameters are restricted to the matrix , which in this particular example, has nine components. Finally, the nine components that would minimize the potential-energy function are obtained. The method will be illustrated using various one-dimensional examples.

### Bars Under Axial Loads

The Rayleigh Ritz method seeks to find an approximate solution to minimize the potential energy of the system. For plane bars under axial loading, the unknown displacement function is . In this case, the potential energy of the system has the following form (See beams under axial loads and energy expressions):

where is Young’s modulus, is the cross sectional area, and is the work done by all the external forces acting on the bar (concentrated forces and distributed loads). To obtain an approximate solution for , a trial function of a finite number of unknowns is first selected. Then, this trial function has to satisfy the “essential” boundary conditions, which are those boundary conditions that would satisfy the physical stability of the system. In the case of bars under axial loads, the trial function has to satisfy the displacement boundary conditions. The following examples show the method for different external loads, and bar configurations.

#### Example 1: A Bar Under Axial Forces

Let a bar with a length and a constant cross sectional area be aligned with the coordinate axis . Assume that the bar is subjected to a horizontal body force (in the direction of the axis ) units of force/unit length where is a known constant. Also, assume that the bar has a constant force of value applied at the end . If the bar is fixed at the end , find the displacement of the bar by directly solving the differential equation of equilibrium. Also, find the displacement function using the Rayleigh Ritz method assuming a polynomial function for the displacement of the degrees (0, 1, 2, 3, and 4). Assume that the bar is linear elastic with Young’s modulus and that the small strain tensor is the appropriate measure of strain. Ignore the effect of Poisson’s ratio.

##### Solution

**Exact Solution**

First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the beams under axial loading section when and are constant is:

Therefore, the exact solution has the form:

where and can be obtained from the boundary conditions:

Therefore, the exact solution for the displacement and the stress are:

View Mathematica Code

Clear[u, x] a = DSolve[{u''[x] == -c*x/(Ee*A), u'[L] == P/(Ee*A), u[0] == 0}, u, x] u = u[x] /. a[[1]] FullSimplify[u] sigma = Ee*D[u, x]; FullSimplify[sigma]

View Python Code

from sympy import * import sympy as sp sp.init_printing(use_latex = "mathjax") x, A, E , C, L, P = symbols("x A E C L P") u = Function("u") u1 = u(x).subs(x,0) u2 = u(x).diff(x).subs(x,L) a = dsolve(u(x).diff(x,2) + C*x/(E*A), u(x), ics = {u1:0, u2:P/(E*A)}) u = a.rhs sigma = E * u.diff(x) display("displacement and stress: ", u, simplify(sigma))

**The Rayleigh Ritz Method:**

The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

Notice that the potential energy lost by the action of the end force is equal to the product of and the displacement evaluated at while the potential energy lost by the action of the distributed body forces is an integral since it acts on each point along the beam length. Further, we will use the constitutive equation to rewrite the potential energy function in terms of the function :

To find an approximate solution, an assumption for the shape or the form of has to be introduced.

**Polynomial of the Zero Degree:**

First, we can assume that is a constant function (a polynomial of the zero degree):

However, to satisfy the “essential” boundary condition that when , the constant has to equal to zero; this leads to the trivial solution , which will automatically be rejected.

**Polynomial of the First Degree:**

As a more refined approximation, we can assume that is a linear function (a polynomial of the first degree):

To satisfy the essential boundary condition @ the coefficient is equal to zero. Therefore:

The strain component can be computed as follows:

Substituting into the equation for the potential energy of the system:

The only variable that can be controlled to minimize the potential energy of the system is . Therefore:

Therefore, the best linear function according to the Rayleigh Ritz method is:

A linear displacement would produce a constant stress:

Unfortunately, this solution is far from being accurate, since it does not satisfy the differential equation of equilibrium (Why?). However, within the only possible or allowed linear shapes, the obtained solution is the minimizer of the potential energy.

**Polynomial of the Second Degree:**

As a more refined approximation, we can assume that is a polynomial of the second degree:

To satisfy the essential boundary condition the coefficient is equal to zero. Therefore:

The strain component can be computed as follows:

Substituting into the equation for the potential energy of the system:

and are the variables that can be controlled to minimize the potential energy of the system. Therefore:

and

Solving the above two equations yields:

Therefore, the best second degree polynomial according to the Rayleigh Ritz method is:

A parabolic displacement would produce a linear stress:

View Mathematica Code

u=a2*x^2+a1*x;

PE=Integrate[(1/2)*Ee*A*D[u,x]^2,{x,0,L}]-P*(u/.x->L)-Integrate[c*x*u,{x,0,L}];

Eq1=D[PE,a2]

Eq2=D[PE,a1]

s=Solve[{Eq1==0,Eq2==0},{a2,a1}]

u=u/.s[[1]]

stress=FullSimplify[Ee*D[u,x]]

View Python Code

from sympy import * import sympy as sp sp.init_printing(use_latex = "mathjax") x, a1, a2, E, A, P, c = symbols("x a_1 a_2 E A P c") u = a2*x**2+a1*x uL = u.subs(x,L) PE = integrate((1/2)*E*A*(u.diff(x)**2), (x,0,L)) - P*uL - integrate(c*x*u, (x,0,L)) display("Potential Energy: ", PE) Eq1 = PE.diff(a2) Eq2 = PE.diff(a1) display("Minimize PE: ", Eq1, Eq2) s = solve((Eq1, Eq2), (a2, a1)) display("Solve: ", s) u = u.subs({a1:s[a1], a2:s[a2]}) display("Best second degree Polynomial(Rayleigh Ritz method): ", u) stress = E * u.diff(x) display("stress: ", simplify(stress))

**Polynomial of the Third Degree:**

Similarly, choosing a polynomial of the third degree as an approximate displacement function means:

with the following approximation for the axial strain

The essential boundary condition lead to . The potential energy of the system has the form:

Clearly, the method generates a large set of equations that are difficult to solve by hand. Using Mathematica we can find the values of the constants , and that would minimize the potential energy of the system:

View Mathematica Code

u=a3*x^3+a2*x^2+a1*x;

PE=Integrate[(1/2)*Ee*A*D[u,x]^2,{x,0,L}]-P*(u/.x->L)-Integrate[c*x*u,{x,0,L}];

Eq1=D[PE,a2]

Eq2=D[PE,a1]

Eq3=D[PE,a3]

s=Solve[{Eq1==0,Eq2==0,Eq3==0},{a2,a1,a3}]

u=u/.s[[1]]

stress=FullSimplify[Ee*D[u,x]]

View Python Code

from sympy import * import sympy as sp sp.init_printing(use_latex = "mathjax") x, a1, a2, a3, E, A, P, c = symbols("x a_1 a_2 a_3 E A P c") u = a3*x**3+a2*x**2+a1*x uL = u.subs(x,L) PE = integrate((1/2)*E*A*(u.diff(x)**2), (x,0,L)) - P*uL - integrate(c*x*u, (x,0,L)) display("Potential Energy: ", PE) Eq1 = PE.diff(a2) Eq2 = PE.diff(a1) Eq3 = PE.diff(a3) display("Minimize PE: ", Eq1, Eq2, Eq3) s = solve((Eq1, Eq2, Eq3), (a2, a1, a3)) display("Solve: ", s) u = u.subs({a1:s[a1], a2:s[a2], a3:s[a3]}) display("Best third degree Polynomial(Rayleigh Ritz method): ", u) stress = E * u.diff(x) display("stress: ", simplify(stress))

Since the assumed trial function is a polynomial of the third degree, which is similar to the exact solution, the Rayleigh Ritz method is able to produce the exact solution for the displacement and the stresses:

#### Example 2:A Bar Under Axial Forces With Varying Cross Sectional Area

Let a bar with a length 2m be aligned with the coordinate axis , and let the width of the bar be equal to 250mm while the height varies linearly such that the height is equal to 500mm when and is equal to 250mm when . Assume that the bar has a constant force of value 200N applied at the end and that the bar is fixed at the end . Find the displacement of the bar by directly solving the differential equation of equilibrium and then using the Rayleigh Ritz method assuming a polynomials of the degrees (0, 1, 2 and 3). Assume the bar to be linear elastic with Young’s modulus and assume that the small strain tensor is the appropriate measure of strain. Ignore the effect of Poisson’s ratio. Compare the solution obtained using the Rayleigh Ritz method with the exact solution.

##### Solution

**Exact Solution:**

First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the beams under axial loading section when is constant and is:

The area varies linearly with according to the equation:

The two boundary conditions are:

Therefore, the exact solution for the displacement and the stress :

Note that the exact solution for the stress could have been directly obtained by dividing the force by the variable cross sectional area. The displacement function would then be obtained by integrating the expression .

View Mathematica Code

Clear[u, x] Ax = 25/100 (5/10 - 25/200*x); Ee = 100000; L = 2; P = 200; s = DSolve[{Ee*u''[x]*Ax + Ee*u'[x]*D[Ax, x] == 0, u'[L] == 200/Ee/(Ax /. x -> L), u[0] == 0}, u, x] u = u[x] /. s[[1]] stress = Ee*D[u, x]

View Python Code

from sympy import * import sympy as sp sp.init_printing(use_latex = "mathjax") x, A, E , L, P = symbols("x A E L P") u = Function("u") u1 = u(x).subs(x,0) u2 = u(x).diff(x).subs(x,L) Ax = 25/100*(5/10-25/200*x) Ax1 = Ax.subs(x,L) display("Area: ", Ax) s = dsolve(E*u(x).diff(x,2)*Ax+E*u(x).diff(x)*Ax.diff(x), u(x), ics = {u1:0, u2:200/E/Ax1}) u = s.rhs.subs({E:100000, L:2, P:200}) stress = (E * u.diff(x)).subs({E:100000, L:2, P:200}) display("displacement and stress: ", u, stress)

**The Rayleigh Ritz Method:**

The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

Notice that the potential energy lost by the action of the end force is equal to the product of and the displacement evaluated at . Further, we will use the constitutive equation to rewrite the potential energy function in terms of the function :

To find an approximate solution, an assumption for the shape or the form of has to be introduced. Similar to the previous example, a polynomial of the zero degree cannot be used as it will automatically render the displacement .

**Polynomial of the First Degree:**

To satisfy the essential boundary condition the coefficient is equal to zero. Therefore:

The strain component can be computed as follows:

Substituting into the equation for the potential energy of the system:

The only variable that can be controlled to minimize the potential energy of the system is . Therefore:

Therefore, the best linear function (units of m) according to the Rayleigh Ritz method is:

A linear displacement would produce a constant stress (units of ):

Unfortunately, this solution is far from being accurate, since it does not satisfy the differential equation of equilibrium (Why?). However, within the only possible or allowed linear shapes, the obtained solution is the minimizer of the potential energy.

**Polynomial of the Second Degree:**

As a more refined approximation, we can assume that is a polynomial of the second degree:

To satisfy the essential boundary condition the coefficient is equal to zero. Therefore:

The strain component can be computed as follows:

Substituting into the equation for the potential energy of the system:

and are the variables that can be controlled to minimize the potential energy of the system. Therefore:

and

Solving the above two equations yields:

Therefore, the best second degree polynomial according to the Rayleigh Ritz method is:

A parabolic displacement would produce a linear stress:

**Polynomial of the Third Degree:**

As a more refined approximation, we can assume that is a polynomial of the third degree:

To satisfy the essential boundary condition the coefficient is equal to zero. Therefore:

The strain component can be computed as follows:

Substituting into the equation for the potential energy of the system:

, and are the variables that can be controlled to minimize the potential energy of the system. Therefore:

Also:

And:

Solving the above three equations yields:

Therefore, the best third degree polynomial according to the Rayleigh Ritz method is:

A cubic displacement would produce a parabolic stress distribution:

View Mathematica Code

Clear[u, u1, u2, u3, u4, stress, stress1, stress2, stress3, x] (*Exact Solution*) Ax = 25/100 (5/10 - 25/200*x); Ee = 100000; L = 2; P = 200; s = DSolve[{Ee*u''[x]*Ax + Ee*u'[x]*D[Ax, x] == 0, u'[L] == 200/Ee/(Ax /. x -> L), u[0] == 0}, u, x] u = u[x] /. s[[1]]; stress = Ee*D[u, x] (*First Degree *) u1 = a1*x; PE = Integrate[1/2*Ee*D[u1, x]^2*Ax, {x, 0, L}] - P*u1 /. x -> L; Eq1 = D[PE, a1]; s = Solve[Eq1 == 0, a1]; u1 = u1 /. s[[1]]; stress1 = Ee*D[u1, x]; (*Second Degree *) u2 = a1*x + a2*x^2; PE = Integrate[1/2*Ee*D[u2, x]^2*Ax, {x, 0, L}] - P*u2 /. x -> L; Eq1 = D[PE, a1]; Eq2 = D[PE, a2]; s = Solve[{Eq1 == 0, Eq2 == 0}, {a1, a2}]; u2 = u2 /. s[[1]]; stress2 = Ee*D[u2, x]; (*Third Degree *) u3 = a1*x + a2*x^2 + a3*x^3; PE = Integrate[1/2*Ee*D[u3, x]^2*Ax, {x, 0, L}] - P*u3 /. x -> L; Eq1 = D[PE, a1]; Eq2 = D[PE, a2]; Eq3 = D[PE, a3]; s = Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0}, {a1, a2, a3}]; u3 = u3 /. s[[1]]; stress3 = Ee*D[u3, x]; (*comparison*) Plot[{u, u1, u2, u3}, {x, 0, L}, PlotLegends -> {Style["Exact", Bold, 12], Style["u1", Bold, 12], Style["u2", Bold, 12], Style["u3", Bold, 12]}, AxesLabel -> {"X1 m", "u(m)"}, BaseStyle -> Directive[Bold, 12]] Plot[{stress, stress1, stress2, stress3}, {x, 0, L}, AxesLabel -> {"X1 m", "stress(Pa)"}, BaseStyle -> Directive[Bold, 12], PlotLegends -> {Style["Exact", Bold, 12], Style["u1", Bold, 12], Style["u2", Bold, 12], Style["u3", Bold, 12]}]

View Python Code

from sympy import * from numpy import * import sympy as sp import numpy as np import matplotlib.pyplot as plt sp.init_printing(use_latex = "mathjax") x, a1, a2, a3, E, P = symbols("x a_1 a_2 a_3 E P") Ax = 25/100*(5/10-25/200*x) display("Exact") u = Function("u") u1_1 = u(x).subs(x,0) u2_1 = u(x).diff(x).subs(x,L) Ax1 = Ax.subs(x,L) display("Area: ", Ax) s = dsolve(E*u(x).diff(x,2)*Ax+E*u(x).diff(x)*Ax.diff(x), u(x), ics = {u1_1:0, u2_1:200/E/Ax1}) u = s.rhs.subs({E:100000, L:2, P:200}) stress0 = (E * u.diff(x)).subs({E:100000, L:2, P:200}) display("displacement and stress: ", u, stress0) display("First Degree") u1 = a1*x uL = u1.subs(x,L) PE = integrate((1/2)*E*((u1.diff(x))**2)*Ax, (x,0,L)) - P*uL PE = PE.subs({E:100000, L:2, P:200}) display("Potential Energy: ", PE) Eq1 = PE.diff(a1) display("Minimize:", Eq1) s = solve(Eq1,a1) display("Solve: ", s) u1 = u1.subs(a1,s[0]) display("Best First degree Polynomial(Rayleigh Ritz method): ", u1) stress1 = (E * u1.diff(x)).subs(E,100000) display("stress: ", stress1) display("Second Degree") u2 = a1*x+a2*x**2 PE2 = integrate((1/2)*E*((u2.diff(x))**2)*Ax, (x,0,L)) - P*u2.subs(x,L) PE2 = PE2.subs({E:100000, L:2, P:200}) display("Potential Energy: ", PE2) Eq1_2 = PE2.diff(a1) Eq2_2 = PE2.diff(a2) s = solve((Eq1_2, Eq2_2),a1, a2) display("Minimize:", Eq1_2, Eq2_2) display("Solve: ", s) u2 = u2.subs({a1:s[a1], a2:s[a2]}) display("Best Second degree Polynomial(Rayleigh Ritz method): ", u2) stress2 = (E * u2.diff(x)).subs(E,100000) display("stress: ", stress2) display("Third Degree") u3 = a1*x+a2*x**2+a3*x**3 PE3 = integrate((1/2)*E*((u3.diff(x))**2)*Ax, (x,0,L)) - P*u3.subs(x,L) PE3 = PE3.subs({E:100000, L:2, P:200}) display("Potential Energy: ", PE3) Eq1_3 = PE3.diff(a1) Eq2_3 = PE3.diff(a2) Eq3_3 = PE3.diff(a3) s = solve((Eq1_3, Eq2_3, Eq3_3),a1, a2, a3) display("Minimize:", Eq1_3, Eq2_3, Eq3_3) display("Solve: ", s) u3 = u3.subs({a1:s[a1], a2:s[a2], a3:s[a3]}) display("Best Third degree Polynomial(Rayleigh Ritz method): ", u3) stress3 = (E * u3.diff(x)).subs(E,100000) display("stress: ", stress3) fig, ax = plt.subplots(1,2, figsize = (15,6)) plt.setp(ax[0], xlabel = "X1 m ", ylabel = "u(m)") plt.setp(ax[1], xlabel = "X1 m", ylabel = "Stress(Pa)") x1 = np.arange(0,2.5,0.5) u_list = [u.subs({x:i}) for i in x1] u1_list = [u1.subs({x:i}) for i in x1] u2_list = [u2.subs({x:i}) for i in x1] u3_list = [u3.subs({x:i}) for i in x1] stress0_list = [stress0.subs({x:i}) for i in x1] stress1_list = [stress1.subs({x:i}) for i in x1] stress2_list = [stress2.subs({x:i}) for i in x1] stress3_list = [stress3.subs({x:i}) for i in x1] display("Comparison") ax[0].plot(x1, u_list, 'blue', label = "exact") ax[0].plot(x1, u1_list, 'orange', label = "u1") ax[0].plot(x1, u2_list, 'green', label = "u2") ax[0].plot(x1, u3_list, 'red', label = "u3") ax[0].legend() ax[1].plot(x1, stress0_list, 'blue', label = "exact") ax[1].plot(x1, stress1_list, 'orange', label = "u1") ax[1].plot(x1, stress2_list, 'green', label = "u2") ax[1].plot(x1, stress3_list, 'red', label = "u3") ax[1].legend()

The exact solution obtained is a logarithmic function, while the approximations used for the Rayleigh Ritz method are polynomial functions. The Taylor series of a logarithmic function can be represented by a linear combination of polynomial functions, and so, the higher the degree of approximation, the closer the approximation to the exact solution. As seen in the plots below, good accuracy is obtained for a polynomial of the second degree for the displacement function; however, the stresses in that case have higher errors. A better approximation for the stresses is obtained using a polynomial of the third degree. In general, the error in the approximate solution (i.e., the error in ) is less than the error in its derivatives which in this case is the value of the strain component . This means as well that in order to achieve higher accuracy in the stress , higher polynomials are needed.

#### Example 3: A Bar Under Varying Body Forces With Fixed Ends

Let a bar with a length 2mbe aligned with the coordinate axis , and let the cross sectional area of the bar be equal to . Assume that the bar is fixed at both ends and . Assume that the bar is subjected to a varying body force that is equal to in the direction of the coordinate axis . Find the displacement of the bar by directly solving the differential equation of equilibrium and then using the Rayleigh Ritz method assuming polynomials of the degrees (0, 1, 2, 3, and 4). Assume that the bar material is linear elastic with Young’s modulus and that the small strain tensor is the appropriate measure of strain. Ignore the effect of Poisson’s ratio. Compare the solution obtained using the Rayleigh Ritz method using three- and four-degree polynomials as approximations with the exact solution.

##### SOLUTION

**Exact Solution:**

Similar to the previous examples, the equilibrium equation is:

Therefore, the solution for has the form:

The boundary conditions for this example are:

Therefore the exact solutions for the displacement (in m) and the stress (in ) are:

**The Rayleigh Ritz Method:**

The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

Notice that the potential energy lost by the action of the distributed body forces is an integral since it acts on each point along the beam length. Further, we will use the constitutive equation to rewrite the potential energy function in terms of the function :

To find an approximate solution, an assumption for the shape or the form of has to be introduced.

**Polynomials of the Zero and First Degrees:**

Assuming that that is a constant function (a polynomial of the zero degree) or a linear function (a polynomial of the first degree) and incorporating the boundary conditions:

leads to the trivial solution , which will automatically be rejected.

**Polynomial of the Second Degree:**

Let be a polynomial of the second degree:

The essential boundary condition can be used to find values for the constants and as follows:

Therefore, the displacement function has only one parameter that can be controlled:

The strain component can be computed as follows:

Substituting into the equation for the potential energy of the system:

can be found by minimizing the potential energy of the system. Therefore:

Therefore, the best second degree polynomial according to the Rayleigh Ritz method is:

A parabolic displacement would produce a linear stress:

**Polynomials of Higher Degrees:**

The results from the following Mathematica code show that the third-degree approximation is very close to the actual solution for both displacements and stresses. The polynomial of the fourth-degree approximation is in fact exact since the exact solution is a polynomial of the fourth degree. Study the Mathematica code carefully to understand how the second boundary conditions of displacement are incorporated. The graph comparing the results are shown below.

View Mathematica Code

p = 5 x^2;

s = DSolve[{Ee*u''[x] == -p/A, u[L] == 0, u[0] == 0}, u, x];

u = u[x] /. s[[1]];

stress = Ee*D[u, x];

(*Second Degree*)

u2 = a2*x^2 + a1*x;

sol = Solve[(u2 /. x -> L) == 0, a2];

u2 = u2 /. sol[[1]];

PE = 1/2*Integrate[Ee*A*D[u2, x]^2, {x, 0, L}] - Integrate[p*u2, {x, 0, L}];

Eq1 = D[PE, a1];

sol1 = Solve[{Eq1 == 0}, {a1}];

u2 = FullSimplify[u2 /. sol1[[1]]];

stress2 = FullSimplify[Ee*D[u2, x]];

(*Third Degree*)

u3 = a3*x^3 + a2*x^2 + a1*x;

sol = Solve[(u3 /. x -> L) == 0, a3];

u3 = u3 /. sol[[1]];

PE = 1/2*Integrate[Ee*A*D[u3, x]^2, {x, 0, L}] - Integrate[p*u3, {x, 0, L}];

Eq1 = D[PE, a1];

Eq2 = D[PE, a2];

sol1 = Solve[{Eq1 == 0, Eq2 == 0}, {a2, a1}];

u3 = FullSimplify[u3 /. sol1[[1]]];

stress3 = FullSimplify[Ee*D[u3, x]];

(*Fourth Degree*)

u4 = a4*x^4 + a3*x^3 + a2*x^2 + a1*x;

sol = Solve[(u4 /. x -> L) == 0, a4];

u4 = u4 /. sol[[1]];

PE = 1/2*Integrate[Ee*A*D[u4, x]^2, {x, 0, L}] - Integrate[p*u4, {x, 0, L}];

Eq1 = D[PE, a1];

Eq2 = D[PE, a2];

Eq3 = D[PE, a3];

sol1 = Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0}, {a1, a2, a3}];

u4 = FullSimplify[u4 /. sol1[[1]]];

stress4 = FullSimplify[Ee*D[u4, x]];

A = (25/100)^2;

Ee = 100000;

L = 2;

(*comparison*)

Plot[{u, u2, u3, u4}, {x, 0, L}, PlotLegends -> {Style["Exact", Bold, 12], Style["u2", Bold, 12], Style["u3", Bold, 12], Style["u4", Bold, 12]}, AxesLabel -> {"X1 m", "u1(m)"}, BaseStyle -> Directive[Bold, 12]]

Plot[{stress, stress2, stress3, stress4}, {x, 0, L}, PlotLegends -> {Style["Exact", Bold, 12], Style["u2", Bold, 12], Style["u3", Bold, 12], Style["u4", Bold, 12]}, AxesLabel -> {"X1 m", "Stress(Pa)"}, BaseStyle -> Directive[Bold, 12]]

View Python Code

from sympy import * from numpy import * import sympy as sp import numpy as np import matplotlib.pyplot as plt sp.init_printing(use_latex = "mathjax") x, a1, a2, a3, a4, E, p, A, L= symbols("x a_1 a_2 a_3 a_4 E p A L") p = 5*x**2 display("Exact") u = Function("u") u1_1 = u(x).subs(x,0) u2_1 = u(x).subs(x,L) s = dsolve(E*u(x).diff(x,2)+ p/A, u(x), ics = {u1_1:0, u2_1:0}) u = s.rhs.subs({E:100000, L:2, A:(25/100)**2}) stress0 = (E * u.diff(x)).subs({E:100000, L:2, A:(25/100)**2}) display("displacement and stress: ", u, stress0) display("Second Degree") u2 = a1*x+a2*x**2 sol = solve(u2.subs(x,L), a2) u2 = u2.subs(a2,sol[0]) display("a2:", sol) display(u2) PE2 = (1/2)*(integrate(E*A*((u2.diff(x))**2), (x,0,L))) - integrate(p*u2, (x,0,L)) PE2 = PE2.subs({E:100000, L:2, A:(25/100)**2}) display("Potential Energy: ", PE2) Eq1_2 = PE2.diff(a1) sol1 = solve(Eq1_2,a1) display("a1:", sol1) u2 = u2.subs(a1,sol1[0]).subs(L,2) display("Best Second degree Polynomial(Rayleigh Ritz method): ", u2) stress2 = (E * u2.diff(x)).subs(E,100000) display("stress: ", stress2) display("Third Degree") u3 = a1*x+a2*x**2+a3*x**3 sol = solve(u3.subs(x,L), a3) display("a3: ", sol) u3 = u3.subs(a3,sol[0]) PE3 = (1/2)*(integrate(E*A*((u3.diff(x))**2), (x,0,L))) - integrate(p*u3, (x,0,L)) PE3 = PE3.subs({E:100000, L:2, A:(25/100)**2}) display("Potential Energy: ", PE3) Eq1_3 = PE3.diff(a1) Eq2_3 = PE3.diff(a2) sol1 = solve((Eq1_3, Eq2_3),a1, a2) display("Minimize:", Eq1_3, Eq2_3) display("Solve: ", sol1) u3 = u3.subs({a1:sol1[a1], a2:sol1[a2], L:2}) display("Best Third degree Polynomial(Rayleigh Ritz method): ", u3) stress3 = (E * u3.diff(x)).subs(E,100000) display("stress: ", stress3) display("Fourth Degree") u4 = a1*x+a2*x**2+a3*x**3+a4*x**4 sol = solve(u4.subs(x,L), a4) display("a4: ", sol) u4 = u4.subs(a4,sol[0]) PE4 = (1/2)*(integrate(E*A*((u4.diff(x))**2), (x,0,L))) - integrate(p*u4, (x,0,L)) PE4 = PE4.subs({E:100000, L:2, A:(25/100)**2}) display("Potential Energy: ", PE4) Eq1_4 = PE4.diff(a1) Eq2_4 = PE4.diff(a2) Eq3_4 = PE4.diff(a3) sol1 = solve((Eq1_4, Eq2_4, Eq3_4),(a1, a2, a3)) display("Minimize:", Eq1_4, Eq2_4, Eq3_4) display("Solve: ", sol1) u4 = u4.subs({a1:sol1[a1], a2:sol1[a2], a3:sol1[a3]}).subs(L,2) display("Best Fourth degree Polynomial(Rayleigh Ritz method): ", u3) stress4 = (E * u4.diff(x)).subs(E,100000) display("stress: ", stress4) fig, ax = plt.subplots(1,2, figsize = (15,6)) plt.setp(ax[0], xlabel = "X1 m ", ylabel = "u(m)") plt.setp(ax[1], xlabel = "X1 m", ylabel = "Stress(Pa)") x1 = np.arange(0,2.1,0.1) u_list = [u.subs({x:i}) for i in x1] u2_list = [u2.subs({x:i}) for i in x1] u3_list = [u3.subs({x:i}) for i in x1] u4_list = [u4.subs({x:i}) for i in x1] stress0_list = [stress0.subs({x:i}) for i in x1] stress2_list = [stress2.subs({x:i}) for i in x1] stress3_list = [stress3.subs({x:i}) for i in x1] stress4_list = [stress4.subs({x:i}) for i in x1] display("Comparison") ax[0].plot(x1, u_list, 'blue', label = "exact") ax[0].plot(x1, u2_list, 'orange', label = "u2") ax[0].plot(x1, u3_list, 'green', label = "u3") ax[0].plot(x1, u4_list, 'red', label = "u4") ax[0].legend() ax[1].plot(x1, stress0_list, 'blue', label = "exact") ax[1].plot(x1, stress2_list, 'orange', label = "u2") ax[1].plot(x1, stress3_list, 'green', label = "u3") ax[1].plot(x1, stress4_list, 'red', label = "u4") ax[1].legend()