パイソニスタの方へ　教えて下さい。
・symbol,varの重複宣言は間違っていますか？
(2024-11-16　ver0.3追加)

・ver0.2 勉強の成果という事で。
・ver0.1
　sympyのPiecewise関数とmaximum関数でチャレンジ。
　??? sympyのPiecewise関数で最大値は計算できませんでした。
　sympyのPiecewise関数で作図plotできました。

オリジナル

基本例題4.16 p95＜「材料力学演習(20221021)一括(ver.3.3)」を勉強したい。#sympy
　　　　　　　　　　？？？タブレット等で、pdfを開く事ができないかも。？？？　

sympyで

・ver0.2
・(未)最大値の計算は、ユーザー定義関数にする予定です。
・【Python3】Decimal.quantizeを用いて小数点第○位以下を四捨五入

# ver0.2
# 基本例題4.16 P95
from sympy import *
from decimal import Decimal, ROUND_HALF_UP

var('FAC,x,a,MAC,FCB,P,MCB,l,RB',real=True)
# var('RA'      ,real=True)
RA = Symbol('RA')
def myRound1(x, d=1):
    p = Decimal(str(x)).quantize(Decimal(str(1/10**d)), rounding=ROUND_HALF_UP)
    p = float(p)
    return p
def myDic (myStrKey,myValue):
    return {sympify(myStrKey):myValue}
FAC=solve(-RA  +FAC        ,FAC)[0]   ;print("# FAC =",FAC)
MAC=solve(-RA*x+MAC        ,MAC)[0]   ;print("# MAC =",MAC)
FCB=solve(-RA  +P+FCB      ,FCB)[0]   ;print("# FAC =",FCB)
MCB=solve(-RA*x+P*(x-a)+MCB,MCB)[0]   ;print("# MCB =",poly(MCB,P).args[0])
RA =solve(MCB.subs({x:l})  ,RA )[0]   ;print("# RA  =",RA)           
RB =solve( RA+RB-P         ,RB )[0]   ;print("# RB  =",RB)                  ;print()
# 
rep = {l:1.0,a:0.4,P:1000}
rep |= myDic("RA",myRound1(RA.subs(rep)))                               
print("# ",rep )
# 
print("# RA  =",myRound1(RA .subs(rep)))
print("# RB  =",myRound1(RB .subs(rep)))
print("# FAC =",myRound1(FAC.subs(rep)))
print("# FCB =",myRound1(FCB.subs(rep)))
print("# MAC =",MAC.subs(rep))
print("# MCB =",MCB.subs(rep))
print("# Fmax =",myRound1(max    (abs(FAC.subs(rep)),abs(FAC.subs(rep)))))
print("# Mmax =",myRound1(maximum(MAC.subs(rep),x,Interval(0.0,0.4)))    )
# FAC = RA
# MAC = RA*x
# FAC = -P + RA
# MCB = P*(a - x) + RA*x
# RA  = P*(-a + l)/l
# RB  = P*a/l

#  {l: 1.0, a: 0.4, P: 1000, RA: 600.0}
# RA  = 600.0
# RB  = 400.0
# FAC = 600.0
# FCB = -400.0
# MAC = 600.0*x
# MCB = 400.0 - 400.0*x
# Fmax = 600.0
# Mmax = 240.0

・ver0.1
　??? エラーです。

# ver0.1
# 基本例題4.16 P95

from sympy import *
l,a,x  =symbols('l a x'  ,real=True)
P,RA,RB=symbols('P RA RB',real=True)
FAC,MAC=symbols('FAC MAC',real=True)
FCB,MCB=symbols('FCB MCB',real=True)
x      =symbols('x'      ,real=True)
#
myF=Eq(-RA+FAC  ,0)
myM=Eq(-RA*x+MAC,0)          
ans=solve([myF,myM],[FAC,MAC])
#
FAC_ans,MAC_ans=ans[FAC],ans[MAC]
myF=Eq(-RA+P+FCB        ,0)
myM=Eq(-RA*x+P*(x-a)+MCB,0)          
ans=solve([myF,myM],[FCB,MCB])
FCB_ans,MCB_ans=ans[FCB],ans[MCB]
RA_ans=solve(MCB_ans.subs({x:l}),RA)[0]
RB_ans=P-RA_ans
#
FAC_ans=FAC_ans.subs({RA:RA_ans})
FCB_ans=FCB_ans.subs({RA:RA_ans})
MAC_ans=MAC_ans.subs({RA:RA_ans})
MCB_ans=MCB_ans.subs({RA:RA_ans})
# print("# RA_ans ",RA_ans )
# print("# RB_ans ",RB_ans )
# print("# FAC_ans",FAC_ans)
# print("# FCB_ans",FCB_ans)
# print("# MAC_ans",MAC_ans)
# print("# MCB_ans",MCB_ans)
# print("################################")
P_v,l_v,a_v=1000,1.0,0.4
mySubs ={P:P_v,l:l_v,a:a_v}
#
print("# RA_ans ",round(RA_ans .subs(mySubs)))
print("# RB_ans ",round(RB_ans .subs(mySubs)))
print("# FAC_ans",round(FAC_ans.subs(mySubs)))
print("# FCB_ans",round(FCB_ans.subs(mySubs)))
print("# MAC_ans",      MAC_ans.subs(mySubs) )
print("# MCB_ans",      MCB_ans.subs(mySubs) )
BMD=Piecewise(( MAC_ans.subs(mySubs) ,x< a),(MCB_ans.subs(mySubs), True)).subs(mySubs)
myMaxBMD=maximum(BMD,x,Interval(0,l_v))  
# plot (BMD,(x,0,l_v))
print("#    BMD ",     BMD)
print("# ????? MaxBMD ",myMaxBMD)
# RA_ans  600
# RB_ans  400
# FAC_ans 600
# FCB_ans -400
# MAC_ans 600.0*x
# MCB_ans 400.0 - 400.0*x
#    BMD  Piecewise((600.0*x, x < 0.4), (400.0 - 400.0*x, True))
#  ????? MaxBMD  0

BMDを作図(sympyのplot moduleで)

(勉強中) Mathematicaで

(勉強中) MATLAB 追加のオプションSymbolic Math Toolbox で

sympyのBeamで

・本命案　勉強中

# ver0.3
# 基本例題4.16 P95
# 作図には、具体的な数値が必要です。
from sympy import *
from sympy.physics.continuum_mechanics.beam import Beam
from sympy import symbols
from decimal import Decimal, ROUND_HALF_UP
E,I,P= symbols('E,I,P')
var("myE,myI,myP,myl,myla")
var("R,x,M,A1,A2,M1,M2,R1,R2")
def remove_last_term(expr):
    # 式を項ごとに分解し、最後の項を除去
    terms = expr.as_ordered_terms()[:-1]  # 最後の項を取り除く
    # 残りの項を使って式を再構成
    return Add(*terms)
def MYshear_stress(expr):
    var('A')
    return remove_last_term(sympify(expr.subs({A:1})))
def myRound1(x, d=1):
    p = Decimal(str(x)).quantize(Decimal(str(1/10**d)), rounding=ROUND_HALF_UP)
    p = float(p)
    return p
def myBeam_Tanjun_reaction_loads(l,E,I,P,la):
    E ,I =symbols('E ,I ')
    b1=Beam(l,E,I)
    E ,I,R1,R2 =symbols('E ,I,R1,R2 ')
    b1.apply_load(R1, 0,-1)
    b1.apply_load(R2, l,-1)
    b1.apply_load(P ,la,-1) 
    # b1.bc_deflection.append([(0,0),(l,0)]) # あってはいけない
    b1.solve_for_reaction_loads(R1,R2)
    R1R2=b1.reaction_loads
    return R1R2[R1],R1R2[R2]
def myBeam_Tanjun_Sakuzu_draw   (l,E,I,P,la):
    E ,I =symbols('E ,I ')
    b2=Beam(l,E,I)
    E ,I,R1,R2 =symbols('E ,I,R1,R2 ')
    b2.apply_load(R1, 0,-1)
    b2.apply_load(R2, l,-1)
    b2.apply_load(P ,la,-1) 
    b2.apply_support(0 ,'pin')
    b2.apply_support(l ,'roller')
    p = b2.draw()  
    p.show() 
    return 
def myBeam_Tanjun_Sakuzu_SFDBMD(l,E,I,P,la):
    E ,I =symbols('E ,I ')
    b3=Beam(l,E,I)
    E ,I,R1,R2 =symbols('E ,I,R1,R2 ')
    b3.apply_load   (R1, 0,-1)
    b3.apply_load   (P ,la,-1) 
    b3.apply_load   (R2, l,-1)
    # b3.bc_deflection.append([(0,0),(l,0)]) # あってはいけない
    # b3.apply_support(0,'pin'   )           # あってはいけない
    # b3.apply_support(l,'roller')           # あってはいけない
    b3.bc_deflection = [(0, 0), (l, 0)]
    b3.solve_for_reaction_loads(R1,R2)
    # 
    print("#",b3.slope())       # ←←←　???　Eの復活
    # b3.plot_slope()
    # plot(b3.slope)
    b3.plot_bending_moment()
    plot(MYshear_stress(b3.shear_stress()),(x,0,l))
    return 
myl,myP,myla=1.0,1000,0.4
R1,R2=myBeam_Tanjun_reaction_loads(myl,myE,myI,myP,myla) ;print("#",myRound1(R1,1),myRound1(R2,1))
myl,myE,myI,myP,myla=1.0,200*(10**9),400*(10**6),1000,0.4
myBeam_Tanjun_Sakuzu_draw  (myl,myE,myI,myP,myla)
myBeam_Tanjun_Sakuzu_SFDBMD(myl,myE,myI,myP,myla)
# -600.0 -400.0
# (-300.0*SingularityFunction(x, 0, 2) + 500*SingularityFunction(x, 0.4, 2) - 200.0*SingularityFunction(x, 1.0, 2) + 64.0)/(E*I)

https://docs.sympy.org/latest/modules/physics/continuum_mechanics/beam.html

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

(テンプレート)

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

sympy.functions.elementary.piecewise.ExprCondPair(expr, cond)

sympy.calculus.util.maximum(f, symbol, domain=Reals)

参考文献

~~＞一部に等分布荷重を受ける単純支持はり~~
~~＞JSME p89~~
　
＞片持ちばり 1)
＞構造力学公式集表5.1 片持ちばりの公式　p128