A well known fact about rocket solid propellant is that it will go off spontaneously provided certain critical conditions of pressure and temperature are met. In order to guarantee  the safety against this possibility a number of verfications are required by the design specifications of its physical layout. One of them addresses the static conditions of the propellant under its own weight and potentially some external load prior to blast-off. The figure presents a model of such configuration. It consists of a variable-section columnar structure loaded with a uniformly distributed load and a concentrated load at its top end. Analyze this structure using analytical/numerical methods and assess the accuracy of your results.

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

1) Define approximate solution (basis functions):

Clear["Global`*"]; 
    Off[General::spell, 
      General::spell1];

Ck[n_] :=  Table[ck[i], {i, n}]
ϕk[n_]  :=  Table[(x^i), {i, n}]
U[n_]  :=  Ck[n] . ϕk[n]
U[1]

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

dUdx[n_]  :=  
      D[U[n], x]  // Simplify
(Subscript[U, tip])[n_]  :=  
      U[n]  /.  x  ->  
          L
q  =  ρg;
A[x_]  :=  
      (Subscript[A, 0]) (2 - (x/( L)))
Π[n_] :=(∫_0^L)( ((EE A[x])/2) (dUdx[n])^2 - 
                  U[n] q) \[DifferentialD]x - P (Subscript[U, tip])[n] // 
        Expand
Π[1]

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

gradΠ[n_]  :=  
      Table[D[Π[n], ck[i]], {i, n}]
gradΠ[1]

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

4) Reconstruct the approximate solution using the obtained coefficients

(Subscript[U, approx])[n_]  :=  U[n]  /.  solution[n]
(Subscript[U, approx])[2]

5) Plot the solution

values  =  {EE  ->  5 (10^9), (Subscript[A, 0]) ->  0.15, 
          L  ->  35,  ρg ->  -900,  
          P ->   -5 (10^3)};
(Subscript[U, exact])  =  (L/(EE (Subscript[A, 0]))) ( ρg x  +  (L ρg  -  P) Log[1 - (x/(2 L))]) /.  values
n  =  2;
(Subscript[U, plot]) =  Evaluate[(Subscript[U, approx])[n]  /.  values]
Plot[{(Subscript[U, exact]), (Subscript[U, plot])}, {x, 0, L /.  values},  
        PlotRange -> All, 
        PlotStyle ->  {{RGBColor[1, 0, 0], Thickness[0.01], 
              Dashing[{0.05, 0.05}]}, Thickness[0.005], 
            Dashing[{0.05, 0.05}]}];

(Subscript[σ, exact]) = EE D[(Subscript[U, exact]), x] /.  values
n = 2;
(Subscript[σ, approx])= EE D[Evaluate[(Subscript[U, approx])[n]], x]/. values

Plot[{(Subscript[σ, exact]), (Subscript[σ, approx])}, {x, 0, L /.  values},  
        PlotRange -> All, 
        PlotStyle ->  {{RGBColor[1, 0, 0], Thickness[0.01], 
              Dashing[{0.05, 0.05}]}, Thickness[0.005], 
            Dashing[{0.05, 0.05}]}];

6) Analyze accuracy of the solution:
Convergence of approximate solutions: tip displacement

δ[
        n_]  :=  (Subscript[U, approx])[n]  /.  x  ->  L /.  values // 
        N
convergence[
        n_]  :=  (δ[n + 1]  -  δ[n])/ δ[
          n + 1]
Table[convergence[i], {i, 4}]

list[n_]   :=  
      Flatten[Table[
          Abs[convergence[i]],  {i, n}]]
ListPlot[Log[10, list[4]], 
        PlotStyle -> PointSize[0.02]];