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



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

:::note info
次　 [(2/3)問題文だけで、図がない場合,AI先生に聞きました。](https://qiita.com/mrrclb48z/items/c325e5459895bc5fe024)
:::

 # オリジナル
(以下は、私が2箇所修正済みです。集中加重→集中荷重。Iy→I)
```
材料力学の問題なのですが・・・
問題：一端固定、他端移動支持の突き出しはりの自由端Bに集中荷重Pが作用している。
このとき、反力RA、RCおよび反モーメントMAを求めよ。また、自由端におけるたわみを求めよ。ただし、
はりの縦弾性係数をE、断面二次モーメントをIとする。
```
画像

https://chie-pctr.c.yimg.jp/dk/iwiz-chie/que-1270054242?w=200&h=200&up=0

https://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1270054242

# 参考文献
- 村上敬宜. (2023). *材料力学* (新装版). 森北出版.　
例題5.13 p109 

# sympyで。Superposition principle で
＞片持ち梁の複数荷重
https://qiita.com/mrrclb48z/items/f2ffdeba7c5d4058c8ea

・点名の付け方が異なります。点名は、知恵袋様です。
・ylからylへの 変換はできませんでした。(応用例5.1参照 P96) 
・l,aの表示と a,bの表示です。
```python
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
```python
# 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行だけです。
```python
# 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 様ありがとうございます。**
・[おまけ(マクスウェルの相反定理) ＜　重ね合わせの原理のすゝめ②~例として京大院の問題を解く~](https://bloodystream.hatenadiary.jp/entry/2021/05/07/193500#:~:text=%E3%81%A8%E6%80%9D%E3%81%84%E3%81%BE%E3%81%99%E3%80%82-,%E3%81%8A%E3%81%BE%E3%81%91(%E3%83%9E%E3%82%AF%E3%82%B9%E3%82%A6%E3%82%A7%E3%83%AB%E3%81%AE%E7%9B%B8%E5%8F%8D%E5%AE%9A%E7%90%86),-%E3%83%9E%E3%82%AF%E3%82%B9%E3%82%A6%E3%82%A7%E3%83%AB%E3%81%AE%E7%9B%B8%E5%8F%8D)
```python
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行だけです。
```python
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で。

![111.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/c6e6dd0c-ac1f-ab6e-ca67-9a1a12af72e2.png)


```python
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の実行環境と 参考のおすすめです。
[(テンプレート)](https://qiita.com/mrrclb48z/items/00dd08b0317069be9342) 　開発環境

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

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

#### wikipedia

https://en.wikipedia.org/wiki/Superposition_principle

https://ja.wikipedia.org/wiki/%E9%87%8D%E3%81%AD%E5%90%88%E3%82%8F%E3%81%9B%E3%81%AE%E5%8E%9F%E7%90%86

#### Qiita

https://qiita.com/mrrclb48z/items/9aa2644f04aa8e46c013