Equivalence of finite element solution using Rayleigh Ritz method and Galerkin method for finding axial deformation of a simple cantilever beam, and on how to determine if the approximate solution found is an exact solution.

by Nasser M. Abbasi  
Student, Engineering Mechanics dept.
University of Wisconsin, Madison
oct 16, 2009

Abstract

Using a simple problem, it is shown that, using a global trial function, the Rayleigh Ritz variational method and the Galerkin method give the same solution. These solutions are compared to the exact solution. It is found that one can determine that the exact solution was reached by increasing the order of the trial function polynomial until the solution returned by Rayleigh Ritz or Galerkin method no longer changes.

The problem

A cantilever beam is loaded with traction load of q = c x (for some constant c). It is fixed on the left side. We need to find axial deformation u (x)

Solving using variational Ritz method

Set up the potential energy (Subscript[π, p]) and Minimize it

({{(Subscript[π, p])=(∫_0^L)(1/2)A E((((du/dx)))^2)dx-(∫_0^L)q u(x) dx }})

where u (x) is the  trial function . The trial function needs to satisfy only the essential conditions, which in this problem  is given as u (0) = 0, so we start by assuming a trial  function u(x).

({{u(x)=(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])(x^2)}})

Since u(0)=0 , we find that (Subscript[a, 0])=0 hence

({{u=(Subscript[a, 1])x+(Subscript[a, 2])(x^2)}})

This trial function has 2 degrees of freedom they are called the generalized  degrees of freedom. The  trial function could have been selected to have only one degree of freedom i.e. u(x)=(Subscript[a, 1])x and this would also have worked but 2 degrees of freedom will give more accurate solution.

Start by defining the trial function

Remove["Global`*"]
u=Sum[(Subscript[a, i]) x^i,{i,0,2}]

Make the trial function satisfy essential boundary condition

sol=First@ Solve[(u/.x->0)==0,(Subscript[a, 0])]

u=u/.sol

Set up the potential energy functional

Π=(∫_0^L)(1/2)A Ε ((((∂_x) u))^2)\[DifferentialD]x-(∫_0^L)c x u \[DifferentialD]x

Solve for ((dΠ )/(Subscript[da, i]))=0for i=1,2 in order to solve for (Subscript[a, 1]) and (Subscript[a, 2])

eq1=(∂_(Subscript[a, 1])) Π==0

eq2=(∂_(Subscript[a, 2])) Π==0

sol=First @Solve[{eq1,eq2},{(Subscript[a, 1]),(Subscript[a, 2])}]

Now that we found (Subscript[a, 1]) and (Subscript[a, 2])we have found our solution

solutionByRitzMethod=u/.%

Solving the problem using Galerkin method

In this method, we start with the differential equation itself and integrate it with a test function we call g (x). The  differential equation for axial deformation with traction load q(x) is  A E(((d^2)u )/(dx^2)) =-q(x), hence our starting point is to solve the integral equation

I=(∫_0^L)(AE(((d^2)u)/(dx^2))+q)g(x)dx=0

We start by assuming a trial function u=(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])(x^2). Now due to u(0)=0 we find that (Subscript[a, 0])=0, hence u=(Subscript[a, 1])x+(Subscript[a, 2])(x^2) Next, we do integration by parts the above integral and obtain

I=[g(x) AE ((Subscript[((du/dx))], 0])^L)-(∫_0^L)(AE(du/dx)+q)((dg(x))/dx)dx+(∫_0^L)q g(x)dx

But

[g(x) AE ((Subscript[((du/dx))], 0])^L)=g(L)AE ((du(L))/dx)-g(0)AE ((du(0))/dx)

and now since a test function must also satisfy the essential conditions, hence g(0)=0 and also the trial function satisfies the natural boundary conditions at x=L which is that there is no load at that point, hence this means A E (du/dx) = 0 since this is the load at the end. Hence the integral becomes

I=-(∫_0^L)(AE(du/dx)+q)((dg(x))/dx)dx+(∫_0^L)q g(x)dx=0

We now pick the test function g(x). In Galerkin method, (Subscript[g, i])(x)=(du/(Subscript[da, i]))hence (Subscript[g, 1])(x)=(du/(Subscript[da, 1]))=x and (Subscript[g, 2])(x)=(du/(Subscript[da, 2]))=(x^2) so now we set up the 2 equations and solve for (Subscript[a, 1]) and (Subscript[a, 2])

u=x (Subscript[a, 1])+(x^2) (Subscript[a, 2]);
(Subscript[g, 1])=x;
(Subscript[g, 2])=(x^2);
q= c x;

eq1=(∫_0^L)(-A Ε (∂_x) u  (∂_x) (Subscript[g, 1])+q (Subscript[g, 1]))\[DifferentialD]x

eq2=(∫_0^L)(-A Ε (∂_x) u  (∂_x) (Subscript[g, 2])+q (Subscript[g, 2]))\[DifferentialD]x

Now solve for (Subscript[a, 1]) and (Subscript[a, 2])

sol=First @Solve[{eq1==0,eq2==0},{(Subscript[a, 1]),(Subscript[a, 2])}]

Substitute the above into the original u(x)

solutionByGalerkinMethod=u/.sol

We see that this solution is the same as the one found from the Ritz method

solutionByGalerkinMethod-solutionByRitzMethod

Exact solution for the problem

The ODE is

({{({{A E u''[x]=-q}})}})

With the essential conditions u(0)=0 and with natural boundary condition A E u'[x]=0 at x=L. Find exact solution:

Remove[u]
q= c x;
sol=First@DSolve[{A Ε u''[x]==-q,u[0]==0,A Ε u'[L]==0},u[x],x];
solutionExact=u[x]/.sol;
Collect[solutionExact,x]

Compare the above solution to the one found above, we see that the approximate solution is

u=((7 c (L^2) x)/(12 A Ε))-((c L (x^2))/(4 A Ε))

while the exact solution is

(Subscript[u, exact=])((c (L^2) x)/(2 A Ε))-((c (x^3))/(6 A Ε))

here is a plot of both solutions. Assume A = 1, c = 1, E = 1 and L = 1, we obtain

Needs["PlotLegends`"]

c=1;Ε=1;L=1;A=1;
Plot[{solutionExact,solutionByRitzMethod},{x,0,L},PlotStyle->{Red,Blue},PlotLegend->{"exact","FEM"},LegendPosition->{-2,-.4}]

How can we improve the FEM solution, can we make it match the exact solution?

Make the trial function a polynomial of higher order. Try

({{u=(Subscript[a, 0])+(Subscript[a, 1])x+(Subscript[a, 2])(x^2)}})+(Subscript[a, 3])(x^3)

and repeat Ritz method. Instead of making the same steps as above, let us write a function which accepts a trial function u(x) and returns back the solution. Find define the function

findSolutionByRitzMethod[u_,x_]:=Module[{eqs,coeffVar,Π},
Clear[A,Ε,c,L];
coeffVar=Flatten[Normal[CoefficientArrays[u,x]]];
coeffVar=Select[coeffVar,UnsameQ[#,0]&];
Π=(∫_0^L)(1/2)A Ε(((D[u,x]))^2)\[DifferentialD]x-(∫_0^L)c x u \[DifferentialD]x;
eqs=Table[ D[Π,coeffVar[[i]] ]==0,{i,1,Length[coeffVar]}];
First @Solve[eqs,coeffVar]
]

Now we can call the above function for the trial function shown above to dtermine the values of the (Subscript[a, i])'s

trialFunction=(Subscript[a, 1])x+(Subscript[a, 2])(x^2)+(Subscript[a, 3])(x^3);
solutionByRitzMethodMoreTerms=findSolutionByRitzMethod[trialFunction,x]

solutionByRitzMethodMoreTerms=trialFunction/.solutionByRitzMethodMoreTerms

We see now that the FEM solution is the same as the exact solution. We can confirm this by taking the difference between them

solutionByRitzMethodMoreTerms-solutionExact//Simplify

What happens if we add more terms to the trial functions? We should expect that the solution would not change. This is confirmed by the following table generated using the above function

numberOfTrialFunctions=5;
(trialFunction=Table[Sum[(Subscript[a, i]) x^i,{i,1,j}],{j,1,numberOfTrialFunctions}])//TableForm

Now find the solution using each of the above trial functions using Ritz method

tbl=Table[findSolutionByRitzMethod[trialFunction[[i]],x],{i,1,numberOfTrialFunctions}];
tbl=Table[{trialFunction[[i]],trialFunction[[i]]/.tbl[[i]]},{i,1,numberOfTrialFunctions}];

Grid[{   {"TrialFunction","Solution"},Sequence@@tbl},Frame->All]

We see from the above that as more terms added to the trial function, after the cubic polynomial, the solution no longer changes. This implies, at least in this example, that the exact solution was reached.  (Need to find a formal proof of this).