The finite element method II

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

Lagrange Polynomials (used as basis functions)

p[n_] := (∏_(i = 1)^n)(x - (Subscript[X, i]))
dp[
        i_, n_]  :=  (D[p[n], x]  /. 
          x -> (Subscript[X, i]))
LagrPoly[i_, n_] :=  
      p[n]/((x - (Subscript[X, i])) dp[i, n])

LagrPoly[1, 3]

ϕ[n_]  := 
      Table[LagrPoly[i, n], {i, n}]
ϕ[3]

Displacement interpolation

UU[n_] :=  Table[U[i], {i, n}]
ue[n_]  := ϕ[n] . UU[n]
ue[2]

Nodal forces

PP[n_] :=  Table[P[i], {i, n}]

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

duedx[n_]  :=  D[ue[n], x]  // Simplify
(Subscript[Π, e])[
        n_]  := (∫_((Subscript[X, 1]))^((Subscript[X, n])))( ((EE A)/2) (duedx[n])^2 - ue[n] (Subscript[q, O])) \[DifferentialD]x  +  PP[n] . UU[n] // 
          Simplify // Expand
(Subscript[Π, e])[2]

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

gradΠ[n_]  :=  
      Table[D[(Subscript[Π, e])[n], U[i]], {i, n}] // 
        FullSimplify
gradΠ[2]

Stiffness[n_]  :=  
      Table[Coefficient[gradΠ[n][[i]], UU[n]], {i, 
          n}]

Stiffness[2]

Det[Stiffness[2]]

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