4610.sagews

# -*- coding: utf-8 -*-
#上面一行表示本程序是采用utf-8字符集编写的。如果去掉，则可能默认是ascii，无法显示汉字

#导入一些packages以及其中的函数
from sympy import *
from sympy.matrices.matrices import *

#定义后面要用到的变量，其中positive=True表示这个变量肯定是正的。如果不显式给出，在开方的时候，可能会出现
#符号i，甚至无法化简(P**2)**(1/2)这种情况
x = Symbol('x')
P = Symbol('P',positive=True)
L = Symbol('L', positive=True)
EI = Symbol('EI', positive=True)

#定义后面要用到的广义坐标,其中 real = True 表示这个变量是实数，理由与上面的类似
a = Symbol('a',real = True)
b = Symbol('b',real = True)
c = Symbol('c',real = True)

#计时用函数
import  time 
START = time.time()

#左边这个井号表示这一行是注释
#下面三行是三种情况对应的命令
#左边是变量得名称，右边依次为： 近似曲线              , [广义位移], 计算方法

y , unknown , METHOD   =   a*(x**2-x**3/(2*L))   , [a]    , "M"

#y , unknown , METHOD  =   a*(x**2-x**3/L)      , [a]      , "y''"
#y ,  unknown , METHOD  =   a*(1-cos(pi*x/2/L)) , [a]      , "M"

#类似的命令只能存在一条，若超过一条，以此框中最后一个为准
#修改此框中命令后，在顶部选择Action...，然后选择evaluate all就会自动计算
#结果见最后两个框

dy = diff(y,x)
dy2 = diff(y,x,2)

#仅仅为了显示结果更好看而已
Stringsub = r"""\begin{array}{l}"""
Stringend = r"""\end{array}"""
Temp = "y =" + latex(y)
Temp += r"\\\frac{\rm{d}y}{\rm{d}x} = "+latex(dy)
Temp += r"\\\frac{{\rm{d}}^2y}{{\rm{d}x}^2} = "+latex(dy2)
String = Stringsub + Temp + Stringend
LatexExpr(String)

\newcommand{\Bold}[1]{\mathbf{#1}}$$\begin{array}{l}y =a \left(x^{2} - \frac{x^{3}}{2 L}\right)\\\frac{\rm{d}y}{\rm{d}x} = a \left(2 x - \frac{3}{2} \frac{x^{2}}{L}\right)\\\frac{{\rm{d}}^2y}{{\rm{d}x}^2} = a \left(2 - 3 \frac{x}{L}\right)\end{array}$$

#对于不同的弹性稳定问题，要修改此框
D = integrate((diff(y,x))**2/2 , (x,0,L))
print "荷载方向位移:Delta = ",D

荷载方向位移:Delta = 17*L**3*a**2/120

#仅仅为了显示结果更好看而已
LatexExpr("\Delta = "+latex(D))

\newcommand{\Bold}[1]{\mathbf{#1}}\Delta = \frac{17}{120} L^{3} a^{2}

if METHOD == "y''":
    PI = integrate(EI*dy2**2/2,(x,0,L))-P*D
elif METHOD == "M":
    M = P*(y.subs(x,L) - y)
    PI = integrate(M**2/2/EI,(x,0,L))-P*D
else:
    print "未定义方法:",METHOD
    PI = 0
print "能量泛函:PI = \n",PI

能量泛函:PI = -17*L**3*P*a**2/120 + 31*L**5*P**2*a**2/(560*EI)

Stringsub = r"\begin{array}{l}"+"\n"
Stringend = "\n"+r"\end{array}"
Temp = r"\Pi = "

if METHOD == "y''":
    Temp += "\n"+r"""\int_0^L 
                        {\frac{1}{2} EI ({ \frac{\rm{d}^2y} {\rm{d}x^2}})^2} 
                  {\rm{d}x}
                   - P\Delta"""
elif METHOD == "M":
    Temp += "\n"+r"""\int_0^L 
                        \frac{1}{2} \frac{M^2}{EI}
                  \rm{d}x
                   - P\Delta"""
else:
    print "未定义方法:",METHOD
    PI = 0
    
Temp += (r"\\"+"\n="+latex(PI))
Temp += Stringend
Temp = Stringsub+Temp
LatexExpr(Temp)

\newcommand{\Bold}[1]{\mathbf{#1}}$$\begin{array}{l} \Pi = \int_0^L \frac{1}{2} \frac{M^2}{EI} \rm{d}x - P\Delta\\ =- \frac{17}{120} L^{3} P a^{2} + \frac{31}{560} \frac{L^{5} P^{2} a^{2}}{EI} \end{array}$$

EQs = Matrix([PI]).jacobian(unknown)
print "平衡方程:"
for e in EQs: print e," = 0",'\n'

平衡方程: -17*L**3*P*a/60 + 31*L**5*P**2*a/(280*EI) = 0

#仅仅为了显示结果更好看而已
Temp = r"""\begin{array}{c}"""
Stringend = "\n"+r"""\end{array}"""

for i in range(len(unknown)):
    Temp += "\n"+r"\\ " + r"\frac{\delta \Pi}{\delta " + latex(unknown[i])+"} =" + latex(EQs[i])+" = 0"
    
Temp += Stringend
Temp = Temp.replace(r"\\","",1)
LatexExpr(Temp)

\newcommand{\Bold}[1]{\mathbf{#1}}$$\begin{array}{c} \frac{\delta \Pi}{\delta a} =- \frac{17}{60} L^{3} P a + \frac{31}{280} \frac{L^{5} P^{2} a}{EI} = 0 \end{array}$$

EQ = Matrix(EQs).jacobian(unknown).det()
print "二阶变分行列式值:\n",pretty(EQ)

二阶变分行列式值: 3 5 2 17⋅L ⋅P 31⋅L ⋅P - ─────── + ──────── 60 280⋅EI

#仅仅为了显示结果更好看而已
LatexExpr(r"|{\delta^2\Pi}| = " + latex(EQ))

\newcommand{\Bold}[1]{\mathbf{#1}}|{\delta^2\Pi}| = - \frac{17}{60} L^{3} P + \frac{31}{280} \frac{L^{5} P^{2}}{EI}

Pcr = solve( EQ, P)
if 0 in Pcr :Pcr.remove(0)
if len(unknown)>1:
    Pcr.sort()
    Model = []
    temp = [solve(EQs.subs(P,min(Pcr)).subs(unknown[0],1)[1:],unknown[1:])]
    for p in Pcr:
        temp = solve(EQs.subs(P,p).subs(unknown[0],1)[1:],unknown[1:])
        m = [1]
        for u in unknown[1:]:
            m.append(temp[u])
        Model.append(m)
else:
    Model = [[1]]
for i in range(len(unknown)):
    print "第",i+1,"模态:",Model[i],"\n对应失稳荷载:Pcr=\n",pretty(Pcr[i]),'\n'

第 1 模态: [1] 对应失稳荷载:Pcr= 238⋅EI ────── 2 93⋅L

#仅仅为了显示结果更好看而已
Stringsub = r"""\begin{array}{l}"""
Stringend = r"""\end{array}"""
Temp = r"P_{cr1} =" + latex(Pcr[0])
for i in range(1,len(Pcr)):
    Temp += r"\\P_{cr" + str(i+1) + "}=" + latex(Pcr[i])
String = Stringsub + Temp + Stringend
LatexExpr(String)

\newcommand{\Bold}[1]{\mathbf{#1}}$$\begin{array}{l}P_{cr1} =\frac{238}{93} \frac{EI}{L^{2}}\end{array}$$

print "数值结果:"
for i in range(len(unknown)):
    print "第",i+1,"模态:",[m.evalf() for m in Model[i]],"\n对应失稳荷载:\n\tPcr=",Pcr[i].evalf(),'\n'

数值结果: 第 1 模态: [1.00000000000000] 对应失稳荷载: Pcr= 2.55913978494624*EI/L**2

print "计算总耗费时间:",time.time()-START,"s"

计算总耗费时间: 6.95691299438 s

"""
    虽然只是一个很简单的题目，手算也不会很难，但是我很重视通用性，希望可以一次性解决所有类似的问题，所以我还是为此写了这段程序，并且调试了一下输出效果。
    有了这段程序，以后就可以集中于位移的计算和能量泛函表达，而不是无聊的推导过程了。也不用一直担心是不是推导过程出错了。因为如果这段程序没有算出结果，问题应该仅仅在于位移函数本身不符合边界条件，或者泛函没写对。
很感谢你的关注。
"""