2.4.2.D.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

two = S(2)
f12 = S(1)/2
f32 = S(3)/2

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

# use constants:
a, q, EI = var('a, q, EI')

# helper factors:
facM = q*a*a
facw = q*a*a*a*a/EI
facwp = q*a*a*a/EI

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

x = var("x")

w1pppp = 1/EI * q
w1ppp = integrate(w1pppp,x) + C1
w1pp = integrate(w1ppp,x) + C2
w1p = integrate(w1pp,x) + C3
w1 =  integrate(w1p,x) + C4

w2pppp = 1/EI * 0
w2ppp = integrate(w2pppp,x) + C5
w2pp = integrate(w2ppp,x) + C6
w2p = integrate(w2pp,x) + C7
w2 =  integrate(w2p,x) + C8

w1 = expand(w1)
w1 = collect(w1, EI)
pprint(latex(w1))
w2 = expand(w2)
w2 = collect(w2, EI)
pprint(latex(w2))


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

# 1: w1(0) = 0
# 2: EI w1''(a) = 0
# 3: w1(a) = w2(0)
# 4: w1'(a) = w2'(0)
# 5: EI w1''(a) = EI w2''(0)
# 6: EI w1'''(a) = EI w2'''(0)
# 7: w2(a) = 0
# 8: w2'(a) = 0

bc1 = Eq(w1.subs(x,0), 0)
bc2 = Eq(EI * w1pp.subs(x,0), 0)
bc3 = Eq(w1.subs(x,a), w2.subs(x,0))
bc4 = Eq(w1p.subs(x,a), w2p.subs(x,0))
bc5 = Eq(EI * w1pp.subs(x,a), EI * w2pp.subs(x,0))
bc6 = Eq(EI * w1ppp.subs(x,a), EI * w2ppp.subs(x,0))
bc7 = Eq(w2.subs(x,a), 0)
bc8 = Eq(w2p.subs(x,a), 0)

eqns = [bc1, bc2, bc3, bc4, bc5, bc6, bc7, bc8]

sol = solve(eqns, unknowns)
sC1 = sol[C1]
sC2 = sol[C2]
sC3 = sol[C3]
sC4 = sol[C4]
sC5 = sol[C5]
sC6 = sol[C6]
sC7 = sol[C7]
sC8 = sol[C8]
for s in [sC1, sC2, sC3, sC4, sC5, sC6, sC7, sC8]:
    pprint(latex(s))

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

w1 = w1.subs({C1: sC1, C2: sC2, C3: sC3, C4: sC4})
w2 = w2.subs({C5: sC5, C6: sC6, C7: sC7, C8: sC8})
for w in [w1, w2]:
    w = expand(w)
    w = collect(w, q/EI)
    pprint(latex(w))


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

w1n = w1/facw
w2n = w2/facw
pprint("\nDimensionslose Absenkungen:")
pprint("\nw / qaaa/EI")
pprint("\nw1:")
pprint(expand(w1n))
pprint("\nw2:")
pprint(expand(w2n))

w1p = diff(w1,x)
w2p = diff(w2,x)
w1pn = w1p/facwp
w2pn = w2p/facwp

# plot=False
print("\n--- e: --------------------------------")
if plot:
    # figure with two subplots:
    f, (ax1, ax2) =  plt.subplots(1, 2, sharey=True)
    f.subplots_adjust(wspace=0)
    tmp = [ax.axhline(ls="-", lw=1, color="black") for ax in [ax1,ax2]]
    # axes-decoration:
    tmp = ax1.set_title("Bereich 1")
    tmp = ax2.set_title("Bereich 2")
    tmp = ax1.set_xlabel("$x_1 / a$")
    tmp = ax2.set_xlabel("$x_2 / a$")
    # xticks:
    tmp = [ax.set_xticks([float(0.25),float(0.5),float(0.75)]) for ax in [ax1,ax2]]

    # Make unitless:
    xi = var("xi")
    [f1,f2] = [X.expand().subs(x/a,xi) for X in [w1n,w2n]]
    [g1,g2] = [X.expand().subs(x/a,xi) for X in [w1pn,w2pn]]

    # Make array for plotting:
    xia = linspace(float(0),float(1))

    [f1,f2] = [fill_array(X,xi,xia) for X in [f1,f2]]
    [g1,g2] = [fill_array(X,xi,xia) for X in [g1,g2]]

    # plot functions:
    tmp = ax1.plot(xia, f1, "b-")
    tmp = ax2.plot(xia, f2, "b-", label="$w \\, / \\, \\frac{qa^4}{EI}$")
    tmp = ax1.plot(xia, g1, "r--")
    tmp = ax2.plot(xia, g2, "r--", label="$w' \\, / \\, \\frac{qa^3}{EI}$")

    # grid on each axis:
    tmp = [ax.grid() for ax in [ax1, ax2]]
    tmp = ax2.legend(loc="best")
    f.subplots_adjust(wspace=0)
    plt.gca().invert_yaxis()
    f.savefig('2.4.2.D_w.svg', transparent=True)

    plt.show()

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

w = w1.subs(x,a)
pprint(latex(w))

E = 210 *giga*pascal
I = 4000 *cm**4
qq = 10**4 *newton/m
aa = 4 *m

w=w.subs({EI: E*I, q:qq, a:aa})

pprint(Num(w/cm,3))

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

w1pp=diff(w1p,x)
w1ppp=diff(w1pp,x)
Q = - EI*w1ppp.subs(x,0)
pprint("Q bei A:")
pprint(Q)
pprint(latex(Q))


w2pp=diff(w2p,x)
w2ppp=diff(w2pp,x)
Q = - EI*w2ppp.subs(x,a)
pprint("Q bei B:")
pprint(Q)
pprint(latex(Q))

--- a: -------------------------------- \frac{C_{1} x^{3}}{6} + \frac{C_{2} x^{2}}{2} + C_{3} x + C_{4} + \frac{q x^{4 }}{24 EI} \frac{C_{5} x^{3}}{6} + \frac{C_{6} x^{2}}{2} + C_{7} x + C_{8} --- b: -------------------------------- - \frac{41 a q}{64 EI} 0 \frac{11 a^{3} q}{96 EI} 0 \frac{23 a q}{64 EI} - \frac{9 a^{2} q}{64 EI} - \frac{5 a^{3} q}{128 EI} \frac{19 a^{4} q}{384 EI} --- c: -------------------------------- \frac{q}{EI} \left(\frac{11 x}{96} a^{3} - \frac{41 a}{384} x^{3} + \frac{x^{4 }}{24}\right) \frac{q}{EI} \left(\frac{19 a^{4}}{384} - \frac{5 x}{128} a^{3} - \frac{9 a^{2 }}{128} x^{2} + \frac{23 a}{384} x^{3}\right) --- d: -------------------------------- Dimensionslose Absenkungen: w / qaaa/EI w1: 3 4 11*x 41*x x ---- - ------ + ----- 96*a 3 4 384*a 24*a w2: 2 3 19 5*x 9*x 23*x --- - ----- - ------ + ------ 384 128*a 2 3 128*a 384*a --- e: --------------------------------

--- f: -------------------------------- \frac{19 a^{4} q}{384 EI} 1.51 --- g: -------------------------------- Q bei A: 41*a*q ------ 64 \frac{41 a}{64} q Q bei B: -23*a*q -------- 64 - \frac{23 a}{64} q