```python
from sympy import *
EI, a, q = var("EI, a, q")
pprint("\nFEM-Solution:")

# 1: Stiffness Matrices:
# Element 1
l = 2*a
l2 = l*l
l3 = l*l*l
K = EI/l3 * Matrix(
    [
    [  4*l2 ,   -6*l ,   2*l2 ,    6*l ,   0   ,   0   ],
    [ -6*l  ,   12   ,  -6*l  ,  -12   ,   0   ,   0   ],
    [  2*l2 ,   -6*l ,   4*l2 ,    6*l ,   0   ,   0   ],
    [  6*l  ,  -12   ,   6*l  ,   12   ,   0   ,   0   ],
    [  0    ,    0   ,   0    ,    0   ,   0   ,   0   ],
    [  0    ,    0   ,   0    ,    0   ,   0   ,   0   ],
    ]
    )
# Element 2
l = a
l2 = l*l
l3 = l*l*l
K += EI/l3 * Matrix(
    [
    [  0   ,   0   ,  0    ,   0   ,  0    ,   0       ],
    [  0   ,   0   ,  0    ,   0   ,  0    ,   0       ],
    [  0   ,   0   ,  4*l2 ,  -6*l ,  2*l2 ,   6*l     ],
    [  0   ,   0   , -6*l  ,  12   , -6*l  , -12       ],
    [  0   ,   0   ,  2*l2 ,  -6*l ,  4*l2 ,   6*l     ],
    [  0   ,   0   ,  6*l  , -12   ,  6*l  ,  12       ],
    ]
    )

# 2: BCs:
p0,w0,p1,w1,p2,w2 = var("ψ₀,w₀,ψ₁,w₁,ψ₂,w₂")
M0,F0,M1,F1,M2,F2 = var("M₀,F₀,M₁,F₁,M₂,F₂")

Mq1, Fq1 = -q/12*a*a, q/2*a
Mq2, Fq2 = -Mq1, Fq1

#              0                         1              2
#                                        qqqqqqqqqqqqqqqq
#              |-------------------------A---------------
u = Matrix([  0,0,                     p1,0,         p2,w2      ] )
f = Matrix([ M0,F0,                    Mq1,Fq1+F1,  Mq2,Fq2     ] )
unks =     [ M0,F0,                    p1,F1,        p2,w2      ]

# 3: Solution:
eq = Eq(K*u, f)
sol = solve(eq, unks)
pprint(sol)

pprint("\nMinimum-Total-Potential-Energy-Principle-Solution:")

# 1: Ansatz:
a0, a1, a2, a3 = var("a0, a1, a2, a3")
b0, b1, b2 = var("b0, b1, b2")

order = 2
if (order == 2):
    # w2 has order 2
    b3, b4 = 0, 0
else:
    # w2 has order 4
    b3, b4 = var("b3, b4")

x1, x2 = var("x1, x2")
w1 = a0 + a1*x1 + a2*x1**2 + a3*x1**3
w1p = diff(w1, x1)
w1pp = diff(w1p, x1)

w2 = b0 + b1*x2 + b2*x2**2 + b3*x2**3 + b4*x2**4
w2p = diff(w2, x2)
w2pp = diff(w2p, x2)

pprint("\nw1 and w2:")
pprint(w1)
pprint(w2)


# 2: Using BCs:
pprint("\nElimination of a0, a1, a2, a3, b0 using BCs:")
# w1(0)=0
e1 = Eq(w1.subs(x1, 0))
# w1'(0)=0
e2 = Eq(w1p.subs(x1, 0))
# w1(2l)=0
e3 = Eq(w1.subs(x1, 2*a))
# w2(0)=0
e4 = Eq(w2.subs(x2, 0))
# w1p(2a)=w2p(0)
e5 = Eq(w1p.subs(x1, 2*a), w2p.subs(x2, 0))

eqns, unks = [e1, e2, e3, e4, e5], [a0, a1, a2, a3, b0]
sol = solve(eqns, unks)
pprint(sol)

sub_list=[
    (a0, sol[a0]),
    (a1, sol[a1]),
    (a2, sol[a2]),
    (a3, sol[a3]),
    (b0, sol[b0]),
    ]

pprint("\nw1 and w2:")
w1 = w1.subs(sub_list)
w2 = w2.subs(sub_list)
pprint(w1)
pprint(w2)

pprint("\nw1'' and w2'':")
w1pp = w1pp.subs(sub_list)
w2pp = w2pp.subs(sub_list)
pprint(w1pp)
pprint(w2pp)

# 3: Using Principle:
pprint("\nU1, U2, Uq:")
i1 = w1pp*w1pp
I1 = integrate(i1, x1)
I1 = I1.subs(x1,2*a) - I1.subs(x1,0)
U1 = EI*I1/2
pprint(U1)

i2 = w2pp*w2pp
I2 = integrate(i2, x2)
I2 = I2.subs(x2,a) - I2.subs(x2,0)
U2 = EI*I2/2
pprint(U2)

i2 = q*w2
I2 = integrate(i2, x2)
I2 = I2.subs(x2,a) - I2.subs(x2,0)
Uq = I2
pprint(Uq)

pprint("\nParameters for U1 + U2 - Uq = Min:")
U = U1 + U2 - Uq

e1 = Eq(diff(U, b1))
e2 = Eq(diff(U, b2))
if (order == 2):
    eqns = [e1, e2]
    unks = [b1, b2]
    sol = solve(eqns, unks)
    sub_list=[
        (b1, sol[b1]),
        (b2, sol[b2]),
        ]
    w2 = w2.subs(sub_list)
else:
    e3 = Eq(diff(U, b3))
    e4 = Eq(diff(U, b4))
    eqns = [e1, e2, e3, e4]
    unks = [b1, b2, b3, b4]
    sol = solve(eqns, unks)
    sub_list=[
        (b1, sol[b1]),
        (b2, sol[b2]),
        (b3, sol[b3]),
        (b4, sol[b4]),
        ]
    w2 = w2.subs(sub_list)

pprint(sol)

pprint("\nw2:")
pprint(w2)

pprint("\nw2(a):")
w2 = w2.subs(x2, a)
pprint(w2)

# FEM-Solution:
# ⎧                               2            4          3             3   ⎫
# ⎪    3⋅a⋅q      -11⋅a⋅q       -a ⋅q       3⋅a ⋅q      -a ⋅q       -5⋅a ⋅q ⎪
# ⎨F₀: ─────, F₁: ────────, M₀: ──────, w₂: ──────, ψ₁: ──────, ψ₂: ────────⎬
# ⎪      8           8            4          8⋅EI        4⋅EI        12⋅EI  ⎪
# ⎩                                                                         ⎭
#
# 2. Minimum-Total-Potential-Energy-Principle-Solution:
#
# w1 and w2:
#                   2        3
# a₀ + a₁⋅x₁ + a₂⋅x₁  + a₃⋅x₁
#                   2
# b₀ + b₁⋅x₂ + b₂⋅x₂
#
# Elimination of a0, a1, a2, a3, b0 using BCs:
# ⎧                  -b₁        b₁        ⎫
# ⎪a₀: 0, a₁: 0, a₂: ────, a₃: ────, b₀: 0⎪
# ⎨                  2⋅a          2       ⎬
# ⎪                            4⋅a        ⎪
# ⎩                                       ⎭
#
# w1 and w2:
#        2        3
#   b₁⋅x₁    b₁⋅x₁
# - ────── + ──────
#    2⋅a         2
#             4⋅a
#              2
# b₁⋅x₂ + b₂⋅x₂
#
# w1'' and w2'':
#   b₁   3⋅b₁⋅x₁
# - ── + ───────
#   a         2
#          2⋅a
# 2⋅b₂
#
# U1, U2, Uq:
#      2
# EI⋅b₁
# ──────
#   a
#          2
# 2⋅EI⋅a⋅b₂
#  3         2
# a ⋅b₂⋅q   a ⋅b₁⋅q
# ─────── + ───────
#    3         2
#
# Parameters for U1 + U2 - Uq = Min:
# ⎧     3          2  ⎫
# ⎪    a ⋅q       a ⋅q⎪
# ⎨b₁: ────, b₂: ─────⎬
# ⎪    4⋅EI      12⋅EI⎪
# ⎩                   ⎭
#
# w2:
#  3         2     2
# a ⋅q⋅x₂   a ⋅q⋅x₂
# ─────── + ────────
#   4⋅EI     12⋅EI
#
# w2(a):
#  4
# a ⋅q
# ────
# 3⋅EI

```

cbd9b7e23a36f5e35313bf0f9cea306bed8d1773
ipynb
Jupyter Notebook
ipynb/TM_2/4_BB/2_BL/2.4.2.K-FEM.ipynb
kassbohm/tm-snippets
5e0621ba2470116e54643b740d1b68b9f28bff12
MIT
ipynb/TM_2/4_BB/2_BL/2.4.2.K-FEM.ipynb
kassbohm/tm-snippets
5e0621ba2470116e54643b740d1b68b9f28bff12
MIT
ipynb/TM_2/4_BB/2_BL/2.4.2.K-FEM.ipynb
kassbohm/tm-snippets
5e0621ba2470116e54643b740d1b68b9f28bff12
MIT
Qwen/Qwen-72B
1. YES
2. YES