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

# 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: --------------------------------")

Ah, Av, Bv, Dh, Dv, Cv = var('Ah, Av, Bv, Dh, Dv, Cv')
unknowns = [Ah, Av, Bv, Dh, Dv, Cv]

qa = q*a
eqns = [
    Eq(0, Ah + Dh),
    Eq(0, Av + Bv - Dv - qa),
    Eq(0, -Bv + two*Dv + f32*qa),
    Eq(0, Dh),
    Eq(0, Dv + Cv - qa),
    Eq(0, f12*qa - Cv),
    ]

sol = solve(eqns, unknowns)
Ah = sol[Ah]
Av = sol[Av]
Bv = sol[Bv]
Cv = sol[Cv]
Dh = sol[Dh]
Dv = sol[Dv]

pprint(latex(Ah))
pprint(latex(Av))
pprint(latex(Bv))
pprint(latex(Cv))
pprint(latex(Dh))
pprint(latex(Dv))

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

x = var('x')

M1 = Av*x
M2 = Av*(a+x) + Bv*x - f12*q*x**2
M3 = Cv*(a-x) - f12*q*(a-x)**2
M1n = M1/facM
M2n = M2/facM
M3n = M3/facM

pprint(expand(M1n))
pprint(expand(M2n))
pprint(expand(M3n))

plot = true
# plot = false

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

    # Make unitless:
    xi = var("xi")
    [f1,f2,f3] = [X.expand().subs(x/a,xi) for X in [M1n,M2n,M3n]]

    # Make array for plotting:
    xia = linspace(float(0),float(1))
    [f1,f2,f3] = [fill_array(X,xi,xia) for X in [f1,f2,f3]]

    # plot functions:
    tmp = ax1.plot(xia, f1, "r--")
    tmp = ax2.plot(xia, f2, "r--")
    tmp = ax3.plot(xia, f3, "r--", label="$M \\, / \\, qa^2$")

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

lhr = M2n.subs(x,a+x).simplify()
rhs = M3n.simplify()
pprint(lhr == rhs)

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


C1, C2, C3, C4, C5, C6 = var('C1, C2, C3, C4, C5, C6')
unknowns = [C1, C2, C3, C4, C5, C6]

w1pp = -M1/EI
w1p = integrate(w1pp,x) + C1
w1 =  integrate(w1p,x) + C2
w2pp = -M2/EI
w2p = integrate(w2pp,x) + C3
w2 =  integrate(w2p,x) + C4
w3pp = -M3/EI
w3p = integrate(w3pp,x) + C5
w3 =  integrate(w3p,x) + C6

for (w,wp) in [(w1p, w1), (w2p, w2), (w3p, w3)]:
    wp = expand(wp)
    wp = collect(wp, q/EI)

    w = expand(w)
    w = collect(w, q/EI)
    pprint(latex(wp))
    pprint(latex(w))


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

# 1: w1(0) = 0
# 2: w1(a) = 0
# 3: w2(0) = 0
# 4: w1'(a) = w2'(0)
# 5: w2(a) = w3(0)
# 6: w3(a) = 0
bc1 = Eq(w1.subs(x,0), 0)
bc2 = Eq(w1.subs(x,a), 0)
bc3 = Eq(w2.subs(x,0), 0)
bc4 = Eq(w1p.subs(x,a), w2p.subs(x,0))
bc5 = Eq(w2.subs(x,a), w3.subs(x,0))
bc6 = Eq(w3.subs(x,a), 0)

sol = solve(eqns, unknowns)
eqns = [bc1, bc2, bc3, bc4, bc5, bc6]

sol = solve(eqns, unknowns)
sC1 = sol[C1]
sC2 = sol[C2]
sC3 = sol[C3]
sC4 = sol[C4]
sC5 = sol[C5]
sC6 = sol[C6]

w1 = w1.subs({C1: sC1, C2: sC2})
w2 = w2.subs({C3: sC3, C4: sC4})
w3 = w3.subs({C5: sC5, C6: sC6})
w1n = w1/facw
w2n = w2/facw
w3n = w3/facw

pprint("\nDimensionslose Absenkungen:")
pprint("\nw / qaaa/EI")
pprint("\nw1:")
pprint(expand(w1n))
pprint("\nw2:")
pprint(expand(w2n))
pprint("\nw3:")
pprint(expand(w3n))

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

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

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

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

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

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

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

    plt.show()

print("\n--- f: --------------------------------")
a1 = w2p.subs(x,a)
a2 = w3p.subs(x,0)
al = abs(a1 - a2)

w = w2.subs(x,a)

# pprint(latex(w))
# pprint(al)

E = 210 *giga*pascal
I = 4000 *cm**4
q_ = 10**4 *newton/m
a_ = 4 *m

w = w.subs({x:a_, EI:E*I, q:q_, a:a_})

pprint(Num(w/cm,3))

al = al.subs({x:a_, EI:E*I, q:q_, a:a_})
ald = degrees(al)
pprint(Num(ald,3))

--- a: -------------------------------- 0 - a q \frac{5 a}{2} q \frac{a q}{2} 0 \frac{a q}{2} --- b: -------------------------------- -x --- a 2 3*x x -1 + --- - ---- 2*a 2 2*a 2 x x --- - ---- 2*a 2 2*a

True --- c: -------------------------------- C_{1} x + C_{2} + \frac{a q x^{3}}{6 EI} C_{1} + \frac{a q x^{2}}{2 EI} C_{3} x + C_{4} + \frac{q}{EI} \left(\frac{a^{2} x^{2}}{2} - \frac{a x^{3}}{4} + \frac{x^{4}}{24}\right) C_{3} + \frac{q}{EI} \left(a^{2} x - \frac{3 a}{4} x^{2} + \frac{x^{3}}{6}\rig ht) C_{5} x + C_{6} + \frac{q}{EI} \left(- \frac{a x^{3}}{12} + \frac{x^{4}}{24}\r ight) C_{5} + \frac{q}{EI} \left(- \frac{a x^{2}}{4} + \frac{x^{3}}{6}\right) --- d: -------------------------------- Dimensionslose Absenkungen: w / qaaa/EI w1: 3 x x - --- + ---- 6*a 3 6*a w2: 2 3 4 x x x x --- + ---- - ---- + ----- 3*a 2 3 4 2*a 4*a 24*a w3: 3 4 5 7*x x x - - ---- - ----- + ----- 8 12*a 3 4 12*a 24*a --- e: --------------------------------

--- f: -------------------------------- 19.0 5.82