forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
11.1.1.1.f.py
90 lines (90 loc) · 3.21 KB
/
11.1.1.1.f.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
from sympy import *
from numpy import *
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
sp.init_printing(use_latex = "mathjax")
x, a1, a2, a3, a4, E, p, A, L= symbols("x a_1 a_2 a_3 a_4 E p A L")
p = 5*x**2
print("Exact")
u = Function("u")
u1_1 = u(x).subs(x,0)
u2_1 = u(x).subs(x,L)
s = dsolve(E*u(x).diff(x,2)+ p/A, u(x), ics = {u1_1:0, u2_1:0})
u = s.rhs.subs({E:100000, L:2, A:(25/100)**2})
stress0 = (E * u.diff(x)).subs({E:100000, L:2, A:(25/100)**2})
print("displacement and stress: ", u, stress0)
print("Second Degree")
u2 = a1*x+a2*x**2
sol = solve(u2.subs(x,L), a2)
u2 = u2.subs(a2,sol[0])
print("a2:", sol)
print(u2)
PE2 = (1/2)*(integrate(E*A*((u2.diff(x))**2), (x,0,L))) - integrate(p*u2, (x,0,L))
PE2 = PE2.subs({E:100000, L:2, A:(25/100)**2})
print("Potential Energy: ", PE2)
Eq1_2 = PE2.diff(a1)
sol1 = solve(Eq1_2,a1)
print("a1:", sol1)
u2 = u2.subs(a1,sol1[0]).subs(L,2)
print("Best Second degree Polynomial(Rayleigh Ritz method): ", u2)
stress2 = (E * u2.diff(x)).subs(E,100000)
print("stress: ", stress2)
print("Third Degree")
u3 = a1*x+a2*x**2+a3*x**3
sol = solve(u3.subs(x,L), a3)
print("a3: ", sol)
u3 = u3.subs(a3,sol[0])
PE3 = (1/2)*(integrate(E*A*((u3.diff(x))**2), (x,0,L))) - integrate(p*u3, (x,0,L))
PE3 = PE3.subs({E:100000, L:2, A:(25/100)**2})
print("Potential Energy: ", PE3)
Eq1_3 = PE3.diff(a1)
Eq2_3 = PE3.diff(a2)
sol1 = solve((Eq1_3, Eq2_3),a1, a2)
print("Minimize:", Eq1_3, Eq2_3)
print("Solve: ", sol1)
u3 = u3.subs({a1:sol1[a1], a2:sol1[a2], L:2})
print("Best Third degree Polynomial(Rayleigh Ritz method): ", u3)
stress3 = (E * u3.diff(x)).subs(E,100000)
print("stress: ", stress3)
print("Fourth Degree")
u4 = a1*x+a2*x**2+a3*x**3+a4*x**4
sol = solve(u4.subs(x,L), a4)
print("a4: ", sol)
u4 = u4.subs(a4,sol[0])
PE4 = (1/2)*(integrate(E*A*((u4.diff(x))**2), (x,0,L))) - integrate(p*u4, (x,0,L))
PE4 = PE4.subs({E:100000, L:2, A:(25/100)**2})
print("Potential Energy: ", PE4)
Eq1_4 = PE4.diff(a1)
Eq2_4 = PE4.diff(a2)
Eq3_4 = PE4.diff(a3)
sol1 = solve((Eq1_4, Eq2_4, Eq3_4),(a1, a2, a3))
print("Minimize:", Eq1_4, Eq2_4, Eq3_4)
print("Solve: ", sol1)
u4 = u4.subs({a1:sol1[a1], a2:sol1[a2], a3:sol1[a3]}).subs(L,2)
print("Best Fourth degree Polynomial(Rayleigh Ritz method): ", u3)
stress4 = (E * u4.diff(x)).subs(E,100000)
print("stress: ", stress4)
fig, ax = plt.subplots(1,2, figsize = (15,6))
plt.setp(ax[0], xlabel = "X1 m ", ylabel = "u(m)")
plt.setp(ax[1], xlabel = "X1 m", ylabel = "Stress(Pa)")
x1 = np.arange(0,2.1,0.1)
u_list = [u.subs({x:i}) for i in x1]
u2_list = [u2.subs({x:i}) for i in x1]
u3_list = [u3.subs({x:i}) for i in x1]
u4_list = [u4.subs({x:i}) for i in x1]
stress0_list = [stress0.subs({x:i}) for i in x1]
stress2_list = [stress2.subs({x:i}) for i in x1]
stress3_list = [stress3.subs({x:i}) for i in x1]
stress4_list = [stress4.subs({x:i}) for i in x1]
print("Comparison")
ax[0].plot(x1, u_list, 'blue', label = "exact")
ax[0].plot(x1, u2_list, 'orange', label = "u2")
ax[0].plot(x1, u3_list, 'green', label = "u3")
ax[0].plot(x1, u4_list, 'red', label = "u4")
ax[0].legend()
ax[1].plot(x1, stress0_list, 'blue', label = "exact")
ax[1].plot(x1, stress2_list, 'orange', label = "u2")
ax[1].plot(x1, stress3_list, 'green', label = "u3")
ax[1].plot(x1, stress4_list, 'red', label = "u4")
ax[1].legend()