12.1.3.1.py

import sympy as sp
from sympy import solve,Matrix,diff,integrate
a0,a1,a2,x,L,EA,p = sp.symbols("a_0 a_1 a_2 x L EA p")
u1,u2,u3 = sp.symbols("u_1 u_2 u_3")
u = a0+a1*x+a2*x**2
print("Quadratic Form of Displacement:")
print("u =", u)
eq1 = u.subs({x:0})
eq2 = u.subs({x:L/2})
eq3 = u.subs({x:L})
s = solve((eq1-u1,eq2-u2,eq3-u3),a0,a1,a2)
print("Solve")
print("a_0 =",s[a0],"a_1 =",s[a1],"a_2 =",s[a2])
u = u.subs(s).expand()
print("Rearranging")
print("u =",u)
N1 = u.coeff(u1)
N2 = u.coeff(u2)
N3 = u.coeff(u3)
print("Shape Functions")
B = Matrix([[diff(N1,x),diff(N2,x),diff(N3,x)]])
print("B =",B.T)
K = EA*Matrix([[integrate(B[i]*B[j], (x,0,L)) for i in range(3)] for j in range(3)])
print("K^e =",K)
Nn = Matrix([[N1,N2,N3]]).T
print("N =",Nn)
fe = Matrix([p*integrate(Nn[i],(x,0,L))for i in range(3)])
print("Nodal Force Vector")
print("f^e =",fe)