Clear["Global`*"];

The Ritz Method

Beam bending using Ritz method
Consider a simply supported beam with a uniformly distributed load of the form: q0. Compute an approximate solution. The properties of the beam are EI, L

Show[Import["beam.jpg", "JPG"]];

1) Seek approximate solutions of the form:

Ck[n_] :=  Table[ck[i], {i, n}]
ϕk[n_] := Table[(x^i)(L - x),{i, n}]
(Subscript[w, 0])[n_]  :=  
      Ck[n] . ϕk[n]

Let's see what these functions look like:

i = 4;
(Subscript[w, 0])[i]
Plot[
        Evaluate[ϕk[i] /.  L ->  1], {x, 0, 1}];

2) Define the approximate potential Π(Ck) for this problem

(Subscript[w, xx])[n_]  :=  (∂_(x, x))(Subscript[w, 0])[n] // 
        Simplify
(Subscript[w, xx])[2]

q[x_]  :=  (Subscript[q, 0])
Π[n_] :=(∫_0^L)( (EI/2) (Subscript[w, xx])[n]^2 - 
                (Subscript[w, 0])[n] q[x]) \[DifferentialD]x // Expand

Π[1]

3) Minimize the potential to find the value of the coefficients that best approximate this problem

gradΠ[n_]  := 
      Table[(∂_(ck[i]))Π[n], {i, 
          n}]
gradΠ[2]

solution[n_]  :=  
      Solve[gradΠ[n]  ==  0 ,  Ck[n]] // 
        Flatten

solution[2]

(Subscript[w, final])[n_]  :=  (Subscript[w, 0])[n]  /.  solution[n]

(Subscript[w, final])[2]

4) Analyze and Plot the solution

values  =  {EI ->  70 (10^9)  ((0.05 (0.3^3))/12), 
          L ->  3, (Subscript[q, 0]) ->  -(10^5)};
Plot[Evaluate[(Subscript[w, final])[3] /.  values], {x, 0, L /.  values},  
        PlotRange -> All, 
        PlotStyle ->  {Thickness[0.01]}];

Beam deflection at x = (L/2), Exact value:

values  =  {EI ->  70 (10^9)  ((0.05 (0.3^3))/12), 
          L ->  3, 
          (Subscript[q, 0]) ->  -(10^5)};
(Subscript[δ, exact]) =  (5/384) (((Subscript[q, 0]) (L^4))/EI) /.  
        values

δ[
        n_]  :=  (Subscript[w, final])[n]  /.  x  ->  (L/2) /.  values // 
        N
convergence[
        n_]  :=  (δ[n + 1]  -  δ[n])/ δ[
          n + 1]
deflections   =  Table[N[δ[i], 15],  {i, 4}]