P1P2片持ち梁の複数荷重(1/4) ひずみエネルギとそれを使った材料力学 基本例題8.33 P289 「材料力学演習(20221021)」をsympyでやってみたい。

:::note info
別解法　[基本例題6.12 P164　はりのたわみと不静定はり 
](https://qiita.com/mrrclb48z/items/7ff5fd4dba41eb0cfdef)
:::

:::note info
次　 [(2/4)Gemini先生に私のコードを解説してもらいました。
](https://qiita.com/mrrclb48z/items/3d9ba6940c91f6e73b27)
:::

# オリジナル
[基本例題8.33 P289 (pdf)](http://zairikiweb.starfree.jp/zai_enshuh/zai_enshuh_ver.3.3.pdf#page=294)＜[「材料力学演習(20221021)一括(ver.3.3)」を勉強したい。#sympy](https://qiita.com/mrrclb48z/items/c3274daf5f7a99cfe027#%E4%BD%9C%E6%A5%AD%E4%B8%AD%E7%A7%81%E3%81%AE%E8%A7%A3%E7%AD%94sympy)
？？？タブレット等で、pdfを開く事ができないかも。？？？　

# sympyで。単位荷重法で。
・(2024-12-16追加)
```python
# p289 例題8.33
from sympy import *
EI,P1,P2,a,b,x =symbols('EI,P1,P2,a,b,x')
def myCollect(stｒ,P1,P2):
    return -collect((sympify(str).expand().simplify()),[P1,P2])    
MAB_str="-a*P1-(a+b)*P2+(P1+P2)*x"                                          ;print("#",MAB_str)
MBC_str=MAB_str                   +"-P1*(x-a)"                              ;print("#",MBC_str)
MAB    =myCollect(MAB_str,P1,P2)                                           #;print("#",MAB_str)
MBC    =myCollect(MBC_str,P1,P2)                                           #;print("#",MBC_str)
δ1     =    1/EI*Integral(-MAB*(-MAB.subs({P1:1,P2:0})),(x,0,a))            ;print("#",δ1)
δ2     =str(1/EI*Integral(-MAB*(-MAB.subs({P1:0,P2:1})),(x,0,a))) + "+" + \
        str(1/EI*Integral(-MBC*(-MBC.subs({P1:0,P2:1})),(x,a,a+b)))         ;print("#",δ2)
# -a*P1-(a+b)*P2+(P1+P2)*x
# -a*P1-(a+b)*P2+(P1+P2)*x-P1*(x-a)
# Integral((-a + x)*(P1*(-a + x) + P2*(-a - b + x)), (x, 0, a))/EI
# Integral((P1*(-a + x) + P2*(-a - b + x))*(-a - b + x), (x, 0, a))/EI+Integral(P2*(-a - b + x)**2, (x, a, a + b))/EI
```
# sympyで。Castigliano's second theoremで。
・出力の順番。注意です。
```python
# ver0.1
from sympy import *
EI,P1,P2,a,b,x =symbols('EI,P1,P2,a,b,x')
MAB_str="-a*P1-(a+b)*P2+(P1+P2)*x"                              #;print("#",MAB_str)
MBC_str=MAB_str                  +"-P1*(x-a)"                   #;print("#",MBC_str)
U      = "1/(2*EI)*integrate((" + MAB_str+")**2,(x,0,a  ))+" \
        +"1/(2*EI)*integrate((" + MBC_str+")**2,(x,a,a+b))"     #;print("#",U)    
print("#",collect((diff(sympify(U),P1).expand()),[P1,P2]))
print("#",collect((diff(sympify(U),P2).expand()),[P1,P2]))
# P2*(a**3/(3*EI) + a**2*b/(2*EI)) + P1*a**3/(3*EI)
# P1*(a**3/(3*EI) + a**2*b/(2*EI)) + P2*(a**3/(3*EI) + a**2*b/EI + a*b**2/EI + b**3/(3*EI))
```
# sympyで。モーメントのつりあいより。
#### 微分方程式の直接解法で

https://qiita.com/mrrclb48z/items/7ff5fd4dba41eb0cfdef#%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%81%AE%E7%9B%B4%E6%8E%A5%E8%A7%A3%E6%B3%95%E3%81%A7

#### The Myosotis Methodで
```python
# ver0.3
from sympy import *
var('EI,P1,P2,a,b')
δ1 =(P2*b)*a**2/(2*EI)+(P1+P2)*a**3/(3*EI)   ;print("#",δ1)
θ1 =(P2*b)*a   /   EI +(P1+P2)*a**2/(2*EI)   ;print("#",θ1) 
δ2d=P2*b**3/(3*EI)                          #;print("#",δ2d)
δ2 =factor(δ1+θ1*b+δ2d)                      ;print("#",δ2)
# P2*a**2*b/(2*EI) + a**3*(P1 + P2)/(3*EI)
# P2*a*b/EI + a**2*(P1 + P2)/(2*EI)
# (2*P1*a**3 + 3*P1*a**2*b + 2*P2*a**3 + 6*P2*a**2*b + 6*P2*a*b**2 + 2*P2*b**3)/(6*EI)
```
#### Castigliano's second theoremで。
```python
# ver0.4
from sympy import *
EI,P1,P2,a,b,x =symbols('EI,P1,P2,a,b,x')
M01_str="-(P1*(a-x)+P2*(a+b-x))"                                ;print("#",M01_str)
M12_str=          "-P2*(a+b-x)"                                 ;print("#",M12_str)
U      = "1/(2*EI)*integrate((" + M01_str+")**2,(x,0,a  ))+" \
        +"1/(2*EI)*integrate((" + M12_str+")**2,(x,a,a+b))"     ;print("#",U)    
δ1=collect((diff(sympify(U),P1).expand()),[P1,P2])              ;print("#",δ1)
δ2=collect((diff(sympify(U),P2).expand()),[P1,P2])              ;print("#",δ2) 
# P2*(a**3/(3*EI) + a**2*b/(2*EI)) + P1*a**3/(3*EI)
# P1*(a**3/(3*EI) + a**2*b/(2*EI)) + P2*(a**3/(3*EI) + a**2*b/EI + a*b**2/EI + b**3/(3*EI))
```
# 割り込みでテスト2例
・以下は単体で、実行できません。ver0.4 の続き
```
# ver0.4 の続き
var('P,L')
rep={P1:P,a:L,P2:0,b:0}
print("#",δ1.subs(rep),diff(δ1.subs(rep),L))
rep={P1:4*P,a:L,P2:P,b:3*L}
print("#",δ1.subs(rep),diff(δ1.subs(rep),L))
print("#",δ2.subs(rep),diff(δ2.subs(rep),L))
# L**3*P/(3*EI) L**2*P/EI
# 19*L**3*P/(6*EI) 19*L**2*P/(2*EI)
# 86*L**3*P/(3*EI) 86*L**2*P/EI
```

# 申し訳ありません。作図は適当にしました。sympyのBeamで。
![111.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/db870618-73dc-e372-dc26-094b5a56aa81.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,P1,P2")
# 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_P1P2_Sakuzu(a,b,E,I,P1,P2):
    b=Beam(a+b,E,I)
    # b.apply_load(PB,a+b,-1)   # error 
    b.apply_load   (P1,0.75,-1) 
    b.apply_load   (P2,1.00,-1) 
    b.apply_support(0,'fixed')
    p = b.draw()  
    p.show() 
    return 
# δ =myBeam_Katamoti_RP_Tawami(a,b,E,I,PB)             
myBeam_Katamoti_P1P2_Sakuzu   (0.75,0.25,1,1,1,1)             
```
# いつもの? sympyの実行環境と 参考のおすすめです。
[(テンプレート)](https://qiita.com/mrrclb48z/items/00dd08b0317069be9342) 　開発環境

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

**「弾性曲線方程式の様々な表示」** の表＜ 各諸量とたわみの関係 ＜ wikipedia

https://ja.wikipedia.org/wiki/%E5%BC%BE%E6%80%A7%E6%9B%B2%E7%B7%9A%E6%96%B9%E7%A8%8B%E5%BC%8F#%E5%90%84%E8%AB%B8%E9%87%8F%E3%81%A8%E3%81%9F%E3%82%8F%E3%81%BF%E3%81%AE%E9%96%A2%E4%BF%82

・上記の英語版は、見つける事ができませんでした。

https://docs.sympy.org/latest/modules/integrals/integrals.html#module-sympy.integrals

#### 参考文献

~~＞一部に等分布荷重を受ける単純支持はり~~
~~＞JSME p89~~
　
~~＞単純支持ばり 7)~~
~~＞構造力学公式集 p136~~