11.1.2.1.a.py

from sympy import *
import sympy as sp 
sp.init_printing(use_latex = "mathjax")
x, EI, P, q, a2, a3, a4 = symbols("x EI P q a_2 a_3 a_4")
print("Exact")
y = Function("y")
th = y(x).diff(x)
M = EI*y(x).diff(x,2)
V = EI*y(x).diff(x,3)
q = 0
M1 = M.subs(x,0)
M2 = M.subs(x,L)
V1 = V.subs(x,0)
V2 = V.subs(x,L)
y1 = y(x).subs(x,0)
y2 = y(x).subs(x,L)
th1 = th.subs(x,0)
th2 = th.subs(x,L)
s = dsolve(EI*y(x).diff(x,4) - q, y(x), ics = {y1:0, th1:0, V2/EI:P/EI, M2/EI:0})
y = s.rhs
th = y.diff(x)
M = EI * y.diff(x,2)
V = EI * y.diff(x,3)
print("y,th,M,V: ", y,th,M,V)
print("Second Degree")
y2 = a2*x**2
PE2 = (EI/2)*integrate(y2.diff(x,2)**2, (x,0,L))+P*y2.subs(x,L)
print("Potential Energy: ", PE2)
Eq2 = PE2.diff(a2)
Sol = solve(Eq2, a2)
print("a2", Sol)
y2 = y2.subs(a2,Sol[0])
th = y2.diff(x)
M = EI * y2.diff(x,2)
V = EI * y2.diff(x,3)
print("y,th,M,V: ", y2,th,M,V)
print("Third Degree")
y3 = a3*x**3+a2*x**2
PE3 = (EI/2)*integrate(y3.diff(x,2)**2, (x,0,L))+P*y3.subs(x,L)
print("Potential Energy: ", PE3)
Eq1_3 = PE3.diff(a3)
Eq2_3 = PE3.diff(a2)
Sol = solve((Eq1_3, Eq2_3), (a2,a3))
print("a2,a3:", Sol)
y3 = y3.subs({a2:Sol[a2], a3:Sol[a3]})
th = y3.diff(x)
M = EI * y3.diff(x,2)
V = EI * y3.diff(x,3)
print("y,th,M,V: ", y3,th,M,V)
print("Fourth Degree")
y4 = a4*x**4+a3*x**3+a2*x**2
PE4 = (EI/2)*integrate(y4.diff(x,2)**2, (x,0,L))+P*y4.subs(x,L)
print("Potential Energy: ", PE4)
Eq1_4 = PE4.diff(a3)
Eq2_4 = PE4.diff(a2)
Eq3_4 = PE4.diff(a4)
Sol = solve((Eq1_4, Eq2_4, Eq3_4), (a2,a3,a4))
print("a2,a3, a4:", Sol)
y4 = y4.subs({a2:Sol[a2], a3:Sol[a3], a4:Sol[a4]})
th = y4.diff(x)
M = EI * y4.diff(x,2)
V = EI * y4.diff(x,3)
print("y,th,M,V: ", y3,th,M,V)