Qiitaで予期しない表示がでています。材料力学「一端固定、他端移動支持の突き出しはりの自由端Bに集中荷重P」(1/3)sympyでやってみたい。

上は予期しない表示です。(204-12-19)
???スライドモードになってしまいました。使い方は、謎でした。

・本ページの「sympyで。Betti's theoremで。」オススメです。(2024-12-19)
・順番がバラバラで申し訳ありません。 (2024-12-19)
・タイトル変更予定。(2024-12-09)
・sympyのcollectを初めて使いました。(2024-12-09)

次　 (2/3)問題文だけで、図がない場合,AI先生に聞きました。

オリジナル

(以下は、私が2箇所修正済みです。集中加重→集中荷重。Iy→I)

材料力学の問題なのですが・・・
問題：一端固定、他端移動支持の突き出しはりの自由端Bに集中荷重Pが作用している。
このとき、反力RA、RCおよび反モーメントMAを求めよ。また、自由端におけるたわみを求めよ。ただし、
はりの縦弾性係数をE、断面二次モーメントをIとする。

画像

参考文献

• 村上敬宜. (2023). 材料力学 (新装版). 森北出版.　
  例題5.13 p109

sympyで。Superposition principle で

＞片持ち梁の複数荷重
https://qiita.com/mrrclb48z/items/f2ffdeba7c5d4058c8ea

・点名の付け方が異なります。点名は、知恵袋様です。
・ylからylへの 変換はできませんでした。(応用例5.1参照 P96)
・l,aの表示と a,bの表示です。

from sympy import *
EI,P,a,b,l,x,RA,RC,MA =symbols('EI,P,a,b,l,x,RA,RC,MA')
var('l3ma')
def myDicSwap(myDic):
    return {list(rep.values())[0]:list(rep.keys())[0]}
y1    = 1/(6*EI)*P *x**2*(3*l-x)
y2    =-1/(6*EI)*RC*x**2*(3*a-x)
y3    =-1/(6*EI)*RC*a**2*(3*x-a)
Rc_sol=solve(Eq(y1         ,-y2).subs({x:a})      ,RC)[0]          #;print("#",RC_sol)
rep   ={3*l-a:l3ma}
RC_sol=Rc_sol.subs(rep)                                             ;print("# RC=",RC_sol.subs(myDicSwap(rep)))
RA_sol=solve(Eq(RA+RC-P    ,0  ).subs({RC:RC_sol}),RA)[0]           ;print("# RA=",RA_sol.subs(myDicSwap(rep)).simplify())
MA_sol=solve(Eq(MA+RC*a-P*l,0  ).subs({RC:RC_sol}),MA)[0]           ;print("# MA=",MA_sol.subs(myDicSwap(rep)))
y1y3  =collect((y1+y3).subs({RC:RC_sol}),l3ma)                     #;print("#",y1y3)
yl    =y1y3.subs({x:l}).simplify().subs(rep).subs(myDicSwap(rep))   ;print("# yl=",yl)
# 
var('my3EI,my2EI')
replab={l:a+b}
RC    =collect((RC_sol.subs(myDicSwap(rep)).subs(replab).expand()),P)
print("# yl=",str(collect((P*a**3)/my3EI+(P*b*a**2)/my2EI,P).subs({P:RC})).replace("my3EI","(3*EI)").replace("my2EI","(2*EI)"),
          "+",str(y1.subs({x:l}).subs(replab)))
# RC= P*(-a + 3*l)/(2*a)
# RA= 3*P*(a - l)/(2*a)
# MA= P*(a - l)/2
# yl= P*(-a*(-a + 3*l)**2 + 4*l**3)/(12*EI)
# yl= P*(1 + 3*b/(2*a))*(a**3/(3*EI) + a**2*b/(2*EI)) + P*(a + b)**3/(3*EI)

sympyで。反力だけなら,Castigliano's second theorem で。

collect() collects common powers of a term in an expression. For example
https://docs.sympy.org/latest/tutorials/intro-tutorial/simplification.html#collect

# ver0.1
from sympy import *
E,I,P =symbols('E,I,P')
var('x,R')
MAC_str="R*(a-x)-P*(a+b-x)"                                   #;print("#",MAC_str)
MCB_str="P*(a+b-x)"                                           #;print("#",MCB_str)
U      = "1/(2*EI)*integrate((" + MAC_str+")**2,(x,0,a  ))+" \
        +"1/(2*EI)*integrate((" + MCB_str+")**2,(x,a,a+b))"   #;print("#",U)    
print("#",collect(solve(Eq(diff(sympify(U),R),0),R)[0],P))
# P*(1 + 3*b/(2*a))

・???無理すれば、計算部分はprint文1行だけです。

# ver0.2
from sympy import *
EI,P,a,b,x,R =symbols('EI,P,a,b,x,R')
print("#",collect(solve(Eq(diff( 1/(2*EI)*integrate((R*(a-x)-P*(a+b-x))**2,(x,0,a)  ) \
                                +1/(2*EI)*integrate((P*(a+b-x)        )**2,(x,a,a+b)) \
                 ,R),0),R)[0],P))
# P*(1 + 3*b/(2*a))

sympyで。Betti's theoremで。

hatena 様ありがとうございます。
・おまけ(マクスウェルの相反定理) ＜　重ね合わせの原理のすゝめ②~例として京大院の問題を解く~

from sympy import *
var('EI,P,a,b,x,RC')
yC1= P *a**3/(3*EI)+P*a**2/(2*EI)*b                          #;print("#",yC1)
yC2=-RC*a**3/(3*EI)                                          #;print("#",yC2)   # ←←←　これがマクスウェルの相反定理
print("#",        solve(Eq(yC1+yC2,0),RC)[0].simplify()   )   
print("#",        solve(Eq(yC1+yC2,0),RC)[0].factor()     )   
print("#",collect(solve(Eq(yC1+yC2,0),RC)[0].factor(),P)  )   
print("#",collect(solve(Eq(yC1+yC2,0),RC)[0]         ,P)  )   
# P + 3*P*b/(2*a)
# P*(2*a + 3*b)/(2*a)
# P*(2*a + 3*b)/(2*a)
# P*(1 + 3*b/(2*a))

・???無理すれば、計算部分はprint文1行だけです。

from sympy import *
var('EI,P,a,b,x,RC')
print("#", collect(solve(Eq(P *a**3/(3*EI)+P*a**2/(2*EI)*b-RC*a**3/(3*EI),0),RC)[0],P) )   
# P*(1 + 3*b/(2*a))

申し訳ありません。作図は適当にしました。sympyのBeamで。

from sympy import *
from sympy.physics.continuum_mechanics.beam import Beam
from sympy import symbols
E,I,RO,RA,PB=symbols('E,I,RO,RA,PB')
a,b         =symbols('a,b',positive=True)
var("x,M,A1,A2,M1,M2")
# def myBeam_Katamoti_RP_Tawami(a,b,E,I,PB):
#     b=Beam(a+b,E,I)
#     b.apply_load( RO,0,  -1)   
#     b.apply_load( RA,a  ,-1) 
#     # b.apply_load(-PB,a+b,-1) 
#     b.apply_load(-PB,2,-1) 
#     b.bc_deflection.append((0,0))
#     b.bc_deflection.append((a,0))
#     b.bc_slope     .append((0,0))
#     b.solve_for_reaction_loads(RO,RA)
#     return b.deflection()
def myBeam_Katamoti_RP_Sakuzu(a,b,E,I,PB):
    b=Beam(a+b,E,I)
    # b.apply_load(PB,a+b,-1)   # error 
    b.apply_load   (PB,1,-1) 
    b.apply_support(0,'fixed')
    b.apply_support(a,'pin')
    p = b.draw()  
    p.show() 
    return 
# δ =myBeam_Katamoti_RP_Tawami(a,b,E,I,PB)             
myBeam_Katamoti_RP_Sakuzu   (0.75,0.25,1,1,1)             

いつもの? sympyの実行環境と 参考のおすすめです。

(テンプレート) 　開発環境

いつもと違うおすすめです。

参考文献

• 村上敬宜. (2023). 材料力学 (新装版). 森北出版.　
  5.8 カスティリアーノの定理 例題5.13 p109

wikipedia

Qiita