2.4.2.E.sagews

from sympy.physics.units import *
from sympy import *
from sympy import N as Num
import matplotlib.pyplot as plt
from numpy import linspace, full_like, array, arange
from mpmath import degrees

# helper function:
def fill_array(f,x,xia):
    if f.is_constant(x):
        res = full_like(xia,f)
    else:
        res = lambdify(x,f,"numpy")
        res = array(res(xia))
    return res

# symbols:
a, MB, EI = var('a, MB, EI')

print("\n--- a: --------------------------------")

# solve for these variables:
C1, C2, C3, C4 = var('C1, C2, C3, C4')
unknowns = [C1, C2, C3, C4]

x = var("x")

wpppp = 1/EI * 0
wppp = integrate(wpppp,x) + C1
wpp = integrate(wppp,x) + C2
wp = integrate(wpp,x) + C3
w =  integrate(wp,x) + C4

for w_ in [wppp, wpp, wp, w]:
    w_ = expand(w_)
    w_ = collect(w_, EI)
    pprint(latex(w_ ))

print("\n--- b: --------------------------------")

# 1: w(0)  = 0
# 2: w'(0) = 0
# 3: w(a)  = 0
# 4: M(a) = MB <=> - EI w''(a) = MB

bc1 = Eq(w.subs(x,0), 0)
bc2 = Eq(wp.subs(x,0), 0)
bc3 = Eq(- EI * wpp.subs(x,a), MB)
bc4 = Eq(w.subs(x,a), 0)

eqns = [bc1, bc2, bc3, bc4]

sol = solve(eqns, unknowns)
sC1 = sol[C1]
sC2 = sol[C2]
sC3 = sol[C3]
sC4 = sol[C4]
for s in [sC1, sC2, sC3, sC4]:
    pprint(latex(s))
print("\n--- c: --------------------------------")

w = w.subs({C1: sC1, C2: sC2, C3: sC3, C4: sC4})
pprint("\nw / ( MB / EI):")
w = expand(w)
w = collect(w, MB/EI)
pprint(latex(w))


print("\n--- d: --------------------------------")

wn = w / ( MB*a*a / EI )
pprint("\nNormalized deflection:")
pprint("\nw / ( MB*a*a / EI ):")
pprint(expand(wn))

wp = diff(w,x)
wpn = wp / ( MB*a / EI )
pprint("\nw' / ( MB*a / EI ):")
pprint(expand(wpn))

xi = var('xi')
f = wn.expand().subs(x/a,xi)
g = wpn.expand().subs(x/a,xi)
# Make array for plotting:
xia = linspace(float(0),float(1))

f = fill_array(f,xi,xia)
g = fill_array(g,xi,xia)

# plot functions:
tmp = plt.plot(xia, f, "b-", label="$w \\, / \\, \\frac{M^{\\mathsf{B}} a^2}{EI}$")
tmp = plt.plot(xia, g, "r--", label="$w' \\, / \\, \\frac{M^{\\mathsf{B}} a}{EI}$")
tmp = plt.grid()
tmp = plt.axhline(0, color='black')
tmp = plt.xlabel("$x / a$")
tmp = plt.legend()
plt.gca().invert_yaxis()

plt.savefig("2.4.2.E_w.svg", transparent=True)
plt.show()

print("\n--- e: --------------------------------")
eq = Eq(wpn, 0)
sol = solve(eq, x)
pprint("\nSolution of quadratic equation:")
pprint(sol)

wpp=diff(wp,x)
pprint(wpp)

print("\n--- f: --------------------------------")

wmax = w.subs(x,2*a/3)

pprint("\nMax. deflection wmax:")
pprint(wmax)

E = 210 *giga*pascal
I = 4000 *cm**4
MB_ = 1000 *kilo*newton*m
a_ = 4 *m

wmax=wmax.subs({EI: E*I, MB:MB_, a:a_})

pprint(Num(wmax/cm,3))

print("\n--- g: --------------------------------")

wppp=diff(wpp,x)

Q = - EI*wppp.subs(x,0)
pprint("\nVertical forces from beam towards ground:")
pprint("at A:")
pprint(Q)

Q = - EI*wppp.subs(x,a)
pprint("at B:")
pprint(-Q)

print("\n--- h: --------------------------------")
M = - EI*wpp.subs(x,0)
pprint(M)

--- a: -------------------------------- C_{1} C_{1} x + C_{2} \frac{C_{1} x^{2}}{2} + C_{2} x + C_{3} \frac{C_{1} x^{3}}{6} + \frac{C_{2} x^{2}}{2} + C_{3} x + C_{4} --- b: -------------------------------- - \frac{3 MB}{2 EI a} \frac{MB}{2 EI} 0 0 --- c: -------------------------------- w / ( MB / EI): \frac{MB}{EI} \left(\frac{x^{2}}{4} - \frac{x^{3}}{4 a}\right) --- d: -------------------------------- Normalized deflection: w / ( MB*a*a / EI ): 2 3 x x ---- - ---- 2 3 4*a 4*a w' / ( MB*a / EI ): 2 x 3*x --- - ---- 2*a 2 4*a

--- e: -------------------------------- Solution of quadratic equation: 2*a [0, ---] 3 /1 3*x\ MB*|- - ---| \2 2*a/ ------------ EI --- f: -------------------------------- Max. deflection wmax: 2 MB*a ----- 27*EI 7.05 --- g: -------------------------------- Vertical forces from beam towards ground: at A: 3*MB ---- 2*a at B: -3*MB ------ 2*a --- h: -------------------------------- -MB ---- 2