(積分)例題1:集中荷重を受ける片持ち梁＜【② 構造力学の「構造」例題】様を参考にsympyで実行した。

Last updated at 2025-01-18Posted at 2025-01-15

パイソニスタの方へ
　アドバイスをいただけると幸いです。
①式「〜＊P/〜」→　文字列「P＊〜/〜」書き換える方法について
　ver0.2 ChatGPT先生に聞いて、無次元変数x/L と Polyで対応しました。 (履歴省略)
　~~ver0.1 自作関数 def myKakeruPw(v,Pw): で対応しました。~~
　~~無次元化があればokかと思いましたが、~~
　~~「５＊L＊P」、「−３＊L＊P」　に対応できていません。→「５＊P＊L」、「−３＊P＊L」~~
~~②無次元化をx/Lで、する方法。~~
③plotの目盛りを記号(PL,L)で表示する方法。
　y軸の数字だけでも非表示できると助かります。

???記号、符号の読替え表が必要かも。
v → w,w → q,L → l,SFDのグラフの表示,BMDのグラフの表示

次　 (積分)例題2:等分布荷重を受ける片持ち梁へ
　　私のsympyの目次へ

オリジナル

Youtube 仮想仕事の原理の教育に関する研究会様　
例題1 片持ち梁1 集中荷重 (2:50〜6:30)
https://youtu.be/NMP0Q78Wl9Q?t=170
＜目次(18:50)
https://youtu.be/NMP0Q78Wl9Q?t=1130
＜フル【② 構造力学の「構造」例題】 (0:00〜19:21)
https://youtu.be/NMP0Q78Wl9Q

sympyのweb上での実行方法

SymPy Live Shellで。
~~FreeCADのマクロは、以下で実行できません。本ページで使っていません。~~
https://qiita.com/mrrclb48z/items/00dd08b0317069be9342#web%E4%B8%8A%E3%81%AE%E5%AE%9F%E8%A1%8C%E6%96%B9%E6%B3%95sympy-live-shell%E3%81%A7

sympyで

・ver0.2
https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.Poly.all_coeffs

# ver0.2
# (積分)例題1:集中荷重を受ける片持ち梁 ＜【② 構造力学の「構造」例題】
from sympy import *
var('x,C1,C2,C3,C4,P,L,EI,xi,x,y')
def myCom_str(first):
    ret= str(Poly(first,L).LC())+'*L**' + str( Poly(first,L).degree() )          #;print(ret)
    return ret
def myMujigennka(v):
    if 'x' in str(v)==True:
       ret=myMujigennka_ari (v)
    else:
       ret=myMujigennka_nasi(v)
    return ret
def myMujigennka_ari(v):
    v_numer =((numer(v).subs({x:L*xi}).expand()))                                #;print(v_numer)
    com     =Poly(v_numer,xi).LC()                                               #;print(com)
    if str(com)[0]=='-':
       com=com*-1                                                                #;print(com)    
    poly=Poly((v_numer/com).expand()).all_coeffs() [::-1]                        #;print(poly) 
    v_str=""
    for i in range(len(poly)):
        if poly[i]!=0:
           if poly[i]==0:
              pass
           elif poly[i]==1:
              v_str=v_str+"+xi**" +str(i)
           else:
              v_str=v_str+"+" + str(poly[i]) + "*xi**" +str(i)           
    v_str=v_str[1:].replace('+-','-').replace('xi','(x/L)')
    v_str=(myCom_str(com)+"/(" + str(denom(v)) +")*("+v_str+")").replace(' ','')   # ;print(v_str)                               
    return v_str
def myMujigennka_nasi(v):
    v_str=(myCom_str(numer(v))+"/(" + str(denom(v))+")").replace(' ','')   # ;print(v_str)                               
    return v_str
EIv4=0                                     #;print(EIv4)
EIv3=integrate(EIv4,x)+     C1             #;print(EIv3)
EIv2=integrate(EIv3,x)+L   *C2             #;print(EIv2)
EIv1=integrate(EIv2,x)+L**2*C3             #;print(EIv1)
EIv =integrate(EIv1,x)+L**3*C4             #;print(EIv )
sol =solve([Eq(1/EI*EIv .subs({x:0}),0),   # 幾何的境界条件
            Eq(1/EI*EIv1.subs({x:0}),0),
            Eq(    -EIv3.subs({x:0}),P),   # 力学的境界条件
            Eq(    -EIv2.subs({x:L}),0),           
           ],[C1,C2,C3,C4])                # ;print(sol)
v   =(1/EI*(EIv.subs({C1:sol[C1],C2:sol[C2],C3:sol[C3],C4:sol[C4]}))).factor()  #;print(w)
print("#",myMujigennka_ari(     v                ))
print("#",myMujigennka((     v   .subs({x:L}))))
print("#",myMujigennka((diff(v,x).subs({x:L}))))
# P*L**3/(6*EI)*(3*(x/L)**2-1*(x/L)**3)
# P*L**3/(3*EI)
# P*L**2/(2*EI)
#
# rep={P:1,L:1,EI:1}
# plot (diff(v,x,3).subs(rep),(x,0,1),ylim=(-1.5,1.5))  # (-)
# plot (diff(v,x,2).subs(rep),(x,0,1))                  # (-)上に凸

・ver0.1 旧いです。

# ver0.1 旧いです。
# (積分)例題1:集中荷重を受ける片持ち梁 ＜【② 構造力学の「構造」例題】
from sympy import *
var('x,C1,C2,C3,C4,P,L,EI')
def myMujigennka(v,com):
    return ( str(com)+"*("
            +str((v/com).expand())                \
                 .replace("x**2/L**2","(x/L)**2") \
                 .replace("x**3/L**3","(x/L)**3") \
                 .replace("x**4/L**4","(x/L)**4") 
            +")").replace(" ","")
def myKakeruPw(v,Pw):
    return Pw+"*"+str(v).replace("*"+Pw,"")
EIv4=0                                     #;print(EIv4)
EIv3=integrate(EIv4,x)+     C1             #;print(EIv3)
EIv2=integrate(EIv3,x)+L   *C2             #;print(EIv2)
EIv1=integrate(EIv2,x)+L**2*C3             #;print(EIv1)
EIv =integrate(EIv1,x)+L**3*C4             #;print(EIv )
sol =solve([Eq(1/EI*EIv .subs({x:0}),0),   # 幾何的境界条件
            Eq(1/EI*EIv1.subs({x:0}),0),
            Eq(    -EIv3.subs({x:0}),P),   # 力学的境界条件
            Eq(    -EIv2.subs({x:L}),0),           
           ],[C1,C2,C3,C4])                # ;print(sol)
v   =(1/EI*(EIv.subs({C1:sol[C1],C2:sol[C2],C3:sol[C3],C4:sol[C4]}))).factor()  #;print(w)
print("#",myKakeruPw(myMujigennka(v,P*L**3/(6*EI)),"P"))                   
print("#",myKakeruPw(     v   .subs({x:L})        ,"P"))
print("#",myKakeruPw(diff(v,x).subs({x:L})        ,"P"))
# P*L**3/(6*EI)*(3*(x/L)**2-(x/L)**3)
# P*L**3/(3*EI)
# P*L**2/(2*EI)
#
rep={P:1,L:1,EI:1}
plot (diff(v,x,3).subs(rep),(x,0,1),ylim=(-1.5,1.5))  # (-)
plot (diff(v,x,2).subs(rep),(x,0,1))                  # (-)上に凸