3.4.E.py.sagews

from sympy.physics.units import *
from sympy import *
half=S(1)/2

print("\n--- User input: ---------------------")

symbolic = True
# symbolic = False

if symbolic:
    # quantity = symbol:
    a = var("a")
    q = var("q")
    F = var("F")
    M = var("M")
    EI = var("EI") # EI = EI
else:
    # quantity =  factor times unit:
    a = 1 *meter
    q = 1 *newton/meter
    F = 1 *newton
    M = 1 *newton*meter
    EI = 1 *newton*meter*meter

# Boundary values:
M0 = -half*q*a*a -M -a*F
Ma = -M

# Variables:
x = var("x")
C1, C2 = var("C1, C2")

# Function M(x) unknown:
Mx = -half*q*x*x + C1*x + C2

# Compute C1, C2:
eq1 = Eq(Mx.subs(x,0),M0)
eq2 = Eq(Mx.subs(x,a),Ma)
sol = solve([eq1, eq2],[C1, C2])
C1s, C2s = sol[C1], sol[C2]

# M(x) known:
Mx = Mx.subs({C1:C1s, C2:C2s})
Mx = Mx.simplify()

pprint("\n\n M(x):")
pprint(Mx)

Mx_alt = -F*(a-x) -M -half*q*(a-x)**2
Mx_alt = Mx_alt.simplify()
pprint("\n\n M(x) (alternative):")
pprint(Mx_alt)

# W = 1 / 2EI \int M*M dx:
W = integrate( half/EI * Mx*Mx, (x, 0, a) )
W = W.simplify()
pprint("\n\n W Stern:")
pprint(W)
# pprint(latex(W))

p = diff(W,M)
w = diff(W,F)

pprint("\n\n phi:")
pprint(p.simplify())

pprint("\n\n w:")
pprint(w.simplify())

--- User input: --------------------- M(x): 2 2 a *q q*x -F*a - M - ---- - ---- + x*(F + a*q) 2 2 M(x) (alternative): 2 q*(a - x) -F*(a - x) - M - ---------- 2 W Stern: / 2 2 3 2 2 4 2\ a*\20*F *a + 60*F*M*a + 15*F*a *q + 60*M + 20*M*a *q + 3*a *q / ----------------------------------------------------------------- 120*EI phi: / 2 \ a*\3*F*a + 6*M + a *q/ ---------------------- 6*EI w: / 2 \ 2 |F*a M a *q| a *|--- + - + ----| \ 3 2 8 / ------------------- EI