1 数值计算用什么作为理工科的社畜，懂计算会计算是一个必不可少的技能，其中尤其是对于土木工程人来说，结构力学、弹塑性力学、计算力学是数值计算中无法逾越的一道坎。由于Matlab简单使用，好学好操作，工科人往往都喜欢使用Matlab来实现数值算法。但是Matlab有几个缺点：
Matlab是非开源的国外商业软件，存在安全问题以及盗版侵权问题；Matlab的安装包极大，对电脑的的要求很高；Matlab的跨平台能力较弱，编写出来的程序往往只能在安装了Matlab的机器上运行。Matlab占用很高为了解决这些缺点，我们可以转而使用python来编写数值计算程序，当前的python版本支持多进程和多线程计算，numpy和sympy等高性能计算模块的开源共享使得python程序的计算性能和速度已经不输于matlab了。且python是开源的免费软件，一直为国内外的大神所维护能够保证性能和安全；python最新版的安装包才40M左右，十分轻量，在低配电脑上也有十分不错的兼容性；python可以在目前已知所有的操作系统上运行，编写一套代码可以在几乎所有的平台上运行。
基于以上的优点，这里强烈推荐一款python的计算模块：Sympy，他能够实现可视化的符号运算，并且与ipython兼容性十分不错，能够输出latex的可视化计算结果，如下图所示。本文将简要介绍Sympy的常用功能，并基于弹性力学给出一个计算模型作为算例，用于演示sympy在理工科的应用实战。
sympy公式可视化呈现2 sympy的安装与使用sympy是一个开源模块，开源地址在github.com/sympy，代码包含详细的功能文档，建议直接fork下载查看。
2.1 安装sympy已经进入了PyPI，可以使用pip或conda直接安装：
代码语言：python复制# Install the sympy modules using pip
pip install sympy
# Install the sympy modules using conda
conda install -c sympy2.2 在jupyter notebook中显示公式ipython的jupyter notebook支持加载mathjax脚本，能够实现可视化展示latex公式。在使用sympy可视化展示公式时，可以直接通过定义符号变量，并进行相关的运算来实现复杂公式的呈现，如下图所示：
简单的latex公式当然也可以直接输出latex代码以嵌入至latex文档：
代码语言：python复制from sympy import *
x, a, b = symbols('x a b')
y = integrate(exp(-x ** 2), (x, a, b))
# 输出latex代码
latex(y)
### output ###
- \frac{\sqrt{\pi} \operatorname{erf}{\left(a \right)}}{2} + \frac{\sqrt{\pi} \operatorname{erf}{\left(b \right)}}{2}3 sympy常用功能3.1 申明变量通过symbols方法将字符串声明为符号变量，。
代码语言：python复制import sympy
# 声明单个变量
x=sympy.symbols('x')
print(x)
# 声明多个变量,以下三个方法都可以
x,y=sympy.symbols(['x','y'])
x,y=sympy.symbols("x,y")
x,y=sympy.symbols("x y")另外在使用symbols申明新的符号变量时，支持latex的上下标语法，如下图所所示：
变量申明的上下标语法3.2 函数表达式（Expr）3.2.1 构造函数代码语言：python复制#### 函数表达式通过变量的运算构造具体函数，或者通过Function函数构造抽象函数。
# 具体函数
f=sympy.sqrt(3*x*y)+x*sympy.sin(y)+y**2+x**3
# 抽象函数
u=sympy.function('u')3.2.2 变量替换和数字赋值代码语言：python复制#### 变量替换与赋值
# expr.subs()可以实现变量替换，替换成数字实现赋值。
g1=f.subs(x,y) # 将f表达式中的x换成y,并将替换的结果赋给g
g2=f.subs({x:2*x,y:2*y}) # 多次替换，字典
g3=f.subs({x:1,y:2})3.2.3 精确求值expr.evalf((n))可以求一个表达式的保留n位有效数字的精确值
代码语言：python复制#### 精确值
# expr.evalf(n)可以求一个表达式的保留n位有效数字的精确值
g3=f.subs({x:1,y:2})
print(g.evalf(4)) # 保留n位有效数字的精确值，8.3593.2.4微分sympy可以实现求微分，方法如下
代码语言：python复制### 微分
# sympy可以实现自动求微分，方法如下
h1=sympy.diff(f,x)
h1=f.diff(x） #同上
h2=sympy.diff(f,x,2,y,1)
# f对x求2次微分，对y求1次微分3.2.5 积分ympy可以实现自动求不定积分和定积分，区别在于是否传入积分上下限
代码语言：python复制#### 积分
# 可以实现自动求不定积分和定积分，区别在于是否传入积分上下限
l1=sympy.integrate(f,x) #不定积分
l2=sympy.integrate(f,(x,1,3)) # 定积分3.3 极限sympy可以实现求极限，注意极限方向
代码语言：python复制##### sympy可以实现求极限，注意极限方向
# 趋于无穷
lim1=sympy.limit(f,x,sympy.oo)
# 趋于0，默认值dir="0",也就是趋于+0
lim2=sympy.limit(f,x,0)
# 趋于0，默认值dir="+"调整为dir="_",也就是趋于-0
lim3=sympy.limit(f,x,0,dir="-")3.4 解方程sympy可以实现解方程，方法是令Expr=0，所以在解方程时，要先构造一个等于0的左端项。返回结果是一个列表，每一项是一个解。如果是方程组，解列表每一项是一个元组，元组对应位置是对应自变量的值。求解方程是要函数是solveset，使用语法是solveset(equation, variable=None, domain=S.Complexes)，分别是等式（默认等于0,），变量，定义域。sp.solveset(E1,x,domain=sp.Reals)
请注意，函数solve也可以用于求解方程式，solve(equations, variables)
代码语言：python复制#### sympy可以实现解方程，方法是令Expr=0，所以在解方程时，要先构造宇哥
#### 等于0的左端项。返回结果是一个列表，每一项是一个解，如果是方程组，解
#### 解列表每一项是一个元组，元组对应位置是对应自变量的值
func=f-3
# 返回f=3时x的值
sympy.solve(func,x)
# x**2+y**2=1,x+y=1
sympy.solve([x**2+y**2-1,x+y-1],[x,y])3.5 泰勒展开（不常见，但要会用）3.5.1 一元展开sympy可以实现泰勒展开，具体函数抽象函数都可以。但是不能对多元函数同时泰勒展开。
代码语言：python复制#### 一元展开
# sympy可以实现泰勒展开，具体函数抽象函数都可以。但是不能对多元函数同时泰勒展开。
# f对x在0处泰勒展开到4阶(把这句话记住，下边四个先后顺序就能记住）
taylor1=sympy.series(f,x,0,4)
# f对x在0处泰勒展开到4阶，去除皮亚诺余项
taylor2=sympy.series(f,x,0,4).remove0
# 抽象函数u对x在0处泰勒展开到4阶
taylor=sympy.series(u(x),x,0,4)3.5.2 多元展开函数的多元泰勒展开可以参考如下的代码。
代码语言：python复制def Taylor_polynomial_sympy(function_expression, variable_list, evaluation_point, degree):
"""
Mathematical formulation reference:
https://math.libretexts.org/Bookshelves/Calculus/Supplemental_Modules_(Calculus)/Multivariable_Calculus/3%3A_Topics_in_Partial_Derivatives/Taylor__Polynomials_of_Functions_of_Two_Variables
:param function_expression: Sympy expression of the function
:param variable_list: list. All variables to be approximated (to be "Taylorized")
:param evaluation_point: list. Coordinates, where the function will be expressed
:param degree: int. Total degree of the Taylor polynomial
:return: Returns a Sympy expression of the Taylor series up to a given degree, of a given multivariate expression, approximated as a multivariate polynomial evaluated at the evaluation_point
"""
from sympy import factorial, Matrix, prod
import itertools
n_var = len(variable_list)
point_coordinates = [(i, j) for i, j in (zip(variable_list, evaluation_point))] # list of tuples with variables and their evaluation_point coordinates, to later perform substitution
deriv_orders = list(itertools.product(range(degree + 1), repeat=n_var)) # list with exponentials of the partial derivatives
deriv_orders = [deriv_orders[i] for i in range(len(deriv_orders)) if sum(deriv_orders[i]) <= degree] # Discarding some higher-order terms
n_terms = len(deriv_orders)
deriv_orders_as_input = [list(sum(list(zip(variable_list, deriv_orders[i])), ())) for i in range(n_terms)] # Individual degree of each partial derivative, of each term
polynomial = 0
for i in range(n_terms):
partial_derivatives_at_point = function_expression.diff(*deriv_orders_as_input[i]).subs(point_coordinates) # e.g. df/(dx*dy**2)
denominator = prod([factorial(j) for j in deriv_orders[i]]) # e.g. (1! * 2!)
distances_powered = prod([(Matrix(variable_list) - Matrix(evaluation_point))[j] ** deriv_orders[i][j] for j in range(n_var)]) # e.g. (x-x0)*(y-y0)**2
polynomial += partial_derivatives_at_point / denominator * distances_powered
return polynomial3.5.3 查看展开项的系数代码语言：python复制taylor.coeff(x) # 查看taylor1中项（x-x0)的系数3.6 e的展开级数并化简代码语言：python复制# e指数函数的级数展开，并化简
f=sp.series(sp.exp(x),x0=1,n=5)
print(f)
### output ###
E + E*(x_a^b - 1) + E*(x_a^b - 1)**2/2 + E*(x_a^b - 1)**3/6 + E*(x_a^b - 1)**4/24 + O((x_a^b - 1)**5, (x_a^b, 1))
sp.expand(f)
# 输出latex代码
sp.latex(sp.expand(f))3.6.1 e的指数级数的展开3.6.2 化简3.7 表达式具体输入值代码语言：python复制# 表达式输入具体值
expr=sp.exp(x)+1
expr3.8 符号化表达式代码语言：python复制# 符号化表达式
str_expr='(x+1)**2'
expr=sp.sympify(str_expr)
expr3.9 极限代码语言：python复制# 极限
sp.Sum(1/x**2,(x,1,sp.oo)).doit()############ 什么意思
sp.limit((1+1/x)**x,x,sp.oo)3.10 计算导数代码语言：python复制# 计算导数
expr=sp.sin(x)
sp.diff(expr,x,2)
# 多元函数偏导
sp.sin(x*y).diff(x,1)3.10.1 对x求两次导3.10.2 对多元函数求偏导3.11 积分运算(integrate)3.11.1 不定积分3.11.2 定积分3.12 解方程代码语言：python复制# 解方程
E1=sp.Eq(x**2+3*x-4,0)
E1
### domain=sp.Reals用于求解方程
# 求解方程是要函数是solveset，
# 使用语法是solveset(equation, variable=None, domain=S.Complexes/Reals #复数集)，
# 分别是等式（默认等于0,），变量，定义域。
sp.solveset(E1,x,domain=sp.Reals)
E2=sp.Eq(x**2+3*x+4,0)
E2
sp.solveset(E2,x,domain=sp.Complexes)
sp.solveset(E2,x,domain=sp.Reals)3.13 解微分方程代码语言：python复制# 解微分方程
# 建立函数变量
f=sp.symbols('y',cls=sp.Function)
E3=sp.Eq(f(x).diff(x)-2*f(x),sp.sin(x))
sp.dsolve(E3,f(x))3.14 矩阵运算代码语言：python复制#### 矩阵运算
# 构造矩阵
sp.Matrix([[1,-1],[2,3],[3,4]])
sp.Matrix([1,2,3])
# 转置
sp.Matrix([1,2,3]).T
A=sp.Matrix([[1,2],[-2,1]])
B=sp.Matrix([[3,4],[-1,2]]).T
A+B
A*B
A**2
B**2 # 得出结论：B转置后B**2结果也会转置EA.tanspose() 为EA的转置矩阵EA.H 为EA的共轭转置矩阵
3.14.1 伴随矩阵代码语言：python复制A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
A
adj_A=A.adjugate()
adj_A4 通过Rayleigh-Ritz法应用板理论计算板的变形一块薄板在两端受到压力时将会出现屈曲现象，板两端受压的力学模型如下图所示。
板的计算理论使用Rayleigh-Ritz法计算板的变形^2，高阶剪切板理论选用Kirchoff板理论，板的挠度表达式如公式所示^1：
\begin{gathered}u_{x}\left( x,y,z\right) =-z\frac{\partial w(x,y)}{\partial x} \\ \begin{gathered}u_{y}\left( x,y,z\right) =-z\frac{\partial w(x,y)}{\partial y} \\ w\left( x,y,z\right) =w(x,y)\end{gathered} \end{gathered} 其中：u_x 和u_y 分别为板单元x方向和y方向的位移，w\left( x,y,z\right) =w(x,y) 表示假定板的挠度沿z方向处的挠度处处相同。
板内部单元的应变为：
\begin{gathered}\begin{gathered}\begin{gathered}\begin{gathered}\epsilon_{xx} =\frac{\partial u_{x}}{\partial x} \\ \epsilon_{yy} =\frac{\partial u_{y}}{\partial y} \end{gathered} \\ \gamma_{xy} =\frac{\partial u_{x}}{\partial y} +\frac{\partial u_{y}}{\partial x} \end{gathered} \\ \gamma_{xz} =\frac{\partial u_{x}}{\partial z} +\frac{\partial u_{z}}{\partial x} \end{gathered} \\ \gamma_{yz} =\frac{\partial u_{y}}{\partial z} +\frac{\partial u_{z}}{\partial y} \end{gathered} 其中，\epsilon_{ii} 表示i 方向上的正应变，\gamma_{ij} 表示i 方向上朝j 方向上的剪应变。板的单元应力为：
\begin{gathered}\begin{gathered}\begin{gathered}\begin{gathered}\begin{gathered}\sigma^{s\pm }_{xx} =\tau^{s} +\left( \lambda^{s} +2\mu^{s} \right) \epsilon^{s\pm }_{xx} +\left( \lambda^{s} +\tau^{s} \right) \epsilon^{s\pm }_{yy} \\ \sigma^{s\pm }_{yy} =\tau^{s} +\left( \lambda^{s} +2\mu^{s} \right) \epsilon^{s\pm }_{yy} +\left( \lambda^{s} +\tau^{s} \right) \epsilon^{s\pm }_{xx} \end{gathered} \\ \sigma^{s\pm }_{xy} =2\left( \mu^{ss} -\tau^{s} \right) \epsilon^{s\pm }_{xy} +\tau^{s} \frac{\partial u^{s\pm }_{x}}{\partial y} \end{gathered} \\ \sigma^{s\pm }_{yx} =2\left( \mu^{ss} -\tau^{s} \right) \epsilon^{s\pm }_{xy} +\tau^{s} \frac{\partial u^{s\pm }_{y}}{\partial x} \end{gathered} \\ \sigma^{s\pm }_{xz} =\tau^{s} \frac{\partial w}{\partial x} \end{gathered} \\ \sigma^{s\pm }_{yz} =\tau^{s} \frac{\partial w}{\partial y} \end{gathered} 其中，\sigma_{ii} 为板中i方向上的正应力，\sigma_{ij} 为板中i方向上朝j方向上的剪应力，\tau^s 为表面剪应力。\mu^{s} 和\mu^{ss} 分别为和板的物理参数有关的超参数。
令：
\begin{gathered}\begin{gathered}\mathbf{\sigma } =\left[ \begin{array}{lllll}\sigma_{xx} &\sigma_{yy} &\sigma_{xy} &\sigma_{xz} &\sigma_{yz} \end{array} \right]^{T} \\ \mathbf{\epsilon } =\left[ \begin{array}{lllll}\epsilon_{xx} &\epsilon_{yy} &\epsilon_{xy} &\epsilon_{xz} &\epsilon_{yz} \end{array} \right]^{T} \end{gathered} \\ \mathbf{\epsilon_{M} } =\left[ \begin{array}{lllll}\frac{\partial^{2} w}{\partial x^{2}} &\frac{\partial^{2} w}{\partial y^{2}} &\frac{\partial^{2} w}{\partial x\partial y} &0&0\end{array} \right] \end{gathered} 则应力公式可以表示为：
\mathbf{\sigma} = \mathbf{B} \cdot \mathbf{\epsilon } + \mathbf{B^s} \cdot \mathbf{\epsilon _M}其中，\mathbf{B} 和\mathbf{B^s} 的表达式见参考文献^1。
构造系统总势能方程：
\delta \left( U-W\right) =0其中，外力做工的表达式和应变能表达式如下所示：
\delta W=\int_{A} (Nxx\frac{\partial w}{\partial x} \frac{\partial \delta w}{\partial x} +N_{yy}\frac{\partial w}{\partial y} \frac{d\delta w}{\partial y} )dA`$\delta U$`表达式其中涉及到所有的变量计算表达式见参考文献^1中公式的16和17，力边界条件见公式18和19。
板的边界条件可以分为3类：SSSS、CCCC、CCSS三种，分别为四端绞支、四端固支和两端绞支和两端固支。具体的表达式见公式21、22和23。使用高阶多项式拟合板的挠度曲线：
\begin{gathered}\begin{gathered}\Phi_{x} (x,y)=\sum^{\infty }_{m=1} \sum^{\infty }_{n=1} \Phi_{xmn} \frac{dX_{m}(x)}{dx} Y_{n}(y),\left( \Phi =\phi ,\theta ,\lambda \right) \\ \Phi_{y} (x,y)=\sum^{\infty }_{m=1} \sum^{\infty }_{n=1} \Phi_{ymn} X_{m}(x)\frac{dY_{n}(y)}{dy} ,\left( \Phi =\phi ,\theta ,\lambda \right) \end{gathered} \\ w(x,y)=\sum^{\infty }_{m=1} \sum^{\infty }_{n=1} W_{mn}X_{m}(x)Y_{n}(y)\end{gathered}其中：\left( \phi_{xmn} ,\theta_{xmn} ,\lambda_{xmn} ,\phi_{ymn} ,\theta_{ymn} ,\lambda_{ymn} ,W_{mn}\right) 为位移形函数的待定系数。X_m(x),Y_n(y) 为位移形函数，应当选为完备函数，如三角函数、多项式函数或小波函数等。在参考文献中，位移形函数选的是三角函数。
根据系统总势能方程，将位移形函数的表达式代入至外力做功和应变能函数中，总势能方程可以表示成如下形式：
\mathbf{K} \cdot \mathbf{\Phi} = 0其中，
\mathbf{\Phi } =[\begin{array}{lllllll}\phi_{xmn} &\theta_{xmn} &\lambda_{xmn} &\phi_{ymn} &\theta_{ymn} &\lambda_{ymn} &W_{mn}\end{array} ]^{T}\mathbf{K} =\left[ \begin{array}{cccc}K_{11}&K_{12}&...&K_{17}\\ K_{21}&K_{22}&...&K_{27}\\ ...&...&...&...\\ K_{71}&K_{72}&...&K_{77}+(e_{1}+e_{2})N_{cr}\end{array} \right] 由于位移形函数为未知量，要使得总势能方程左侧等于0，则矩阵\mathbf{K} 必须不满秩，即矩阵K的行列式等于0，即\det \left( \mathbf{K} \right)=0
基于以上的推导过程，编写相应的计算程序求解：
代码语言：python复制import numpy as np
from sympy import *
# Double integral
def integrate2(f, x, y):
g = integrate(f, x)
g = integrate(g, y)
return g
# Double differential
def diffn(f, x, n):
while n != 0:
f = diff(f, x)
n = n - 1
return f
# The program is used to solve the problem of critical buckling force of thin plates
ST = 0
# Define the condition of the boundary
BC = 3
m = n = 1
C11 = 263E9
C12 = 154E9
C44 = 127E9
h = 0.8
a = 10
b = 10
# h = 50E-9
# b = a = 10 * h
E = 25.5E9
# E = (C11 - C12) * (C11 + 2 * C12) / (C11 + C12)
v = 0.2
# v = C12 / (C11 + C12)
# Modified coefficient of the shear stress
k = Symbol('k')
x, y, z = symbols('x y z')
w_xx, w_yy, w_xy, w_x, w_y = symbols('w_xx w_yy w_xy w_x w_y')
theta_xx, theta_yy, theta_xy, theta_yx, theta_x, theta_y = symbols('theta_xx theta_yy theta_xy theta_yx theta_x theta_y')
lambda_xx, lambda_yy, lambda_xy, lambda_yx, lambda_x, lambda_y = symbols('lambda_xx lambda_yy lambda_xy lambda_yx lambda_x lambda_y')
phi_xx, phi_yy, phi_xy, phi_yx, phi_x, phi_y = symbols('phi_xx phi_yy phi_xy phi_yx phi_x phi_y')
Gxy = G = E / 2 / (1 + v)
# Gxy = G = C44
Diff = E * h ** 3 / 12 / (1 - v ** 2)
DF = E * h ** 3
NCr = k * pi ** 2 * Diff / (a ** 2)
kv = 0
print('Plate Thickness: h=', h, 'm')
print('Length to thickness ratio: a/h=', a / h)
print('length to width ratio: a/b=', a / b)
if BC == 1:
print('Boundary Condition: SSSS')
# First-order buckling in two directions
alphaM = m * pi / a
beltaN = n * pi / b
Xm = sin(alphaM * x)
Yn = sin(beltaN * y)
D1Xm = cos(x * alphaM) * alphaM ** 1
D2Xm = -sin(x * alphaM) * alphaM ** 2
D3Xm = -cos(x * alphaM) * alphaM ** 3
D4Xm = sin(x * alphaM) * alphaM ** 4
D1Yn = cos(y * beltaN) * beltaN ** 1
D2Yn = -sin(y * beltaN) * beltaN ** 2
D3Yn = -cos(y * beltaN) * beltaN ** 3
D4Yn = sin(y * beltaN) * beltaN ** 4
if BC == 2:
print('Boundary Condition: CCCC')
alphaM = (m + 0.5) * pi / a
beltaN = (n + 0.5) * pi / b
Xm = sin(alphaM * x) - sinh(alphaM * x) - (sin(alphaM * a) - sinh(alphaM * a)) /\
(cos(alphaM * a) - cosh(alphaM * a)) * (cos(alphaM * x) - cosh(alphaM * x))
Yn = sin(beltaN * y) - sinh(beltaN * y) - (sin(beltaN * b) - sinh(beltaN * b)) /\
(cos(beltaN * b) - cosh(beltaN * b)) * (cos(beltaN * y) - cosh(beltaN * y))
D1Xm = cos(x * alphaM) * alphaM - cosh(x * alphaM) * alphaM - (sin(a * alphaM) - sinh(a * alphaM))\
* (-sin(x * alphaM) - sinh(x * alphaM) * alphaM) / (cos(a * alphaM) - cosh(a * alphaM))
D2Xm = -sin(x * alphaM) * alphaM ** 2 - sinh(x * alphaM) * alphaM ** 2 - (sin(a * alphaM) - sinh(a * alphaM)) \
* (-cos(x * alphaM) * alphaM ** 2 - cosh(x * alphaM) * alphaM ** 2) / (cos(a * alphaM) - cosh(a * alphaM))
D3Xm = -cos(x * alphaM) * alphaM **3 - cosh(x * alphaM) * alphaM ** 3 - (sin(a * alphaM) - sinh(a * alphaM)) \
* (sin(x * alphaM) * alphaM ** 3 - sinh(x * alphaM) * alphaM ** 3) / (cos(a * alphaM) - cosh(a * alphaM))
D4Xm = sin(x * alphaM) * alphaM ** 4 - sinh(x * alphaM) * alphaM ** 4 - (sin(a * alphaM) - sinh(a * alphaM)) \
* (cos(x * alphaM) * alphaM ** 4 - cosh(x * alphaM) * alphaM ** 4) / (cos(a * alphaM) - cosh(a * alphaM))
D1Yn = cos(y * beltaN) * beltaN - cosh(y * beltaN) * beltaN - (sin(b * beltaN) - sinh(b * beltaN))\
* (-sin(y * beltaN) * beltaN - sinh(y * beltaN) * beltaN) / (cos(b * beltaN) - cosh(b * beltaN))
D2Yn = -sin(y * beltaN) * beltaN ** 2 - sinh(y * beltaN) * beltaN ** 2 - (sin(b * beltaN) - sinh(b * beltaN)) \
* (-cos(y * beltaN) * beltaN ** 2 - cosh(y * beltaN) * beltaN ** 2) / (cos(b * beltaN) - cosh(b * beltaN))
D3Yn = -cos(y * beltaN) * beltaN ** 3 - cosh(y * beltaN) * beltaN ** 3 - (sin(b * beltaN) - sinh(b * beltaN)) \
* (sin(y * beltaN) * beltaN ** 3 - sinh(y * beltaN) * beltaN ** 3) / (cos(b * beltaN) - cosh(b * beltaN))
D4Yn = sin(y * beltaN) * beltaN ** 4 - sinh(y * beltaN) * beltaN ** 4 - (sin(b * beltaN) - sinh(b * beltaN)) \
* (cos(y * beltaN) * beltaN ** 4 - cosh(y * beltaN) * beltaN ** 4) / (cos(b * beltaN) - cosh(b * beltaN))
if BC == 3:
print('Boundary Conditions: CCSS')
alphaM = (m + 0.5) * pi / a
beltaN = n * pi / b
Xm = sin(alphaM * x) - sinh(alphaM * a) - (sin(alphaM * a) - sinh(alphaM * a)) / (cos(alphaM * a) - cosh(alphaM * a))\
* (cos(alphaM * x) - cosh(alphaM * x))
Ym = sin(beltaN * y)
D1Xm = cos(x * alphaM) * alphaM - cosh(x * alphaM) * alphaM - (sin(a * alphaM) - sinh(a * alphaM))\
* ((-sin(x * alphaM) * alphaM - sinh(x * alphaM) * alphaM)) / (cos(a * alphaM) - cosh(a * alphaM))
D2Xm = -sin(x * alphaM) * alphaM ** 2 - sinh(x * alphaM) * alphaM ** 2 - (sin(a * alphaM) - sinh(a * alphaM))\
* (-cos(x * alphaM) * alphaM ** 2 - cosh(x * alphaM) * alphaM ** 2) / (cos(a * alphaM) - cosh(a * alphaM))
D3Xm = -cos(x * alphaM) * alphaM ** 3 - cosh(x * alphaM) * alphaM ** 3 - (sin(a * alphaM) - sinh(a * alphaM))\
* (sin(x * alphaM) * alphaM ** 3 - sinh(a * alphaM) * alphaM ** 3) / (cos(a * alphaM) - cosh(a * alphaM))
D4Xm = sin(x * alphaM) * alphaM ** 4 - sinh(x * alphaM) * alphaM ** 4 - (sin(a * alphaM) - sinh(a * alphaM))\
* (cos(x * alphaM) * alphaM ** 4 - cosh(x * alphaM) * alphaM ** 4) / (cos(a * alphaM) - cosh(a * alphaM))
D1Yn = cos(y * beltaN) * beltaN
D2Yn = -sin(y * beltaN) * beltaN ** 2
D3Yn = -cos(y * beltaN) * beltaN ** 3
D4Yn = sin(y * beltaN) * beltaN ** 4
PT = 1
for i in range(0, 1, 2):
kar = 1
if PT == 2:
kar = 5 / 6
R1 = R2 = R3 = 0
Rz1 = Rz2 = Rz3 = 0
Rp1 = Rp2 = Rp3 = 0
Rn1 = Rn2 = Rn3 = 0
Rz1p = Rz2p = Rz3p = 0
Rz1n = Rz2n = Rz3n = 0
print('Plate Theory: Kirchoff Plate')
########################## strain ###################################
varepsilon_xx = (R1 - z) * w_xx + R1 * phi_xx + R2 * theta_xx + R3 * lambda_xx
varepsilon_yy = (R1 - z) * w_yy + R1 * phi_yy + R2 * theta_yy + R3 * lambda_yy
gamma_xy = 2 * (R1 - z) * w_xy + R1 * (phi_xy + phi_yx) + R2 * (theta_xy + theta_yx) + R3 * (lambda_xy + lambda_yx)
gamma_xz = Rz1 * (w_x + phi_x) + Rz2 * theta_x + Rz3 * lambda_x
gamma_yz = Rz1 * (w_y + phi_y) + Rz2 * theta_y + Rz3 * lambda_y
########################## corrected stress ###################################
sigma_xx = E / (1 - v ** 2) * (varepsilon_xx + v * varepsilon_yy)
sigma_yy = E / (1 - v ** 2) * (varepsilon_yy + v * varepsilon_xx)
sigma_xy = G * gamma_xy
sigma_xz = G * gamma_xz
sigma_yz = G * gamma_yz
########################## Strain on the top surface ###################################
varepsilon_xxsp = (Rp1 - h / 2) * w_xx + Rp1 * phi_xx + Rp2 * theta_xx + Rp3 * lambda_xx
varepsilon_yysp = (Rp1 - h / 2) * w_yy + Rp1 * phi_yy + Rp2 * theta_yy + Rp3 * lambda_yy
gamma_xysp = 2 * (Rp1 - h / 2) * w_xy + Rp1 * (phi_xy + phi_yx) + Rp2 * (theta_xy + theta_yx) + Rp3 * (lambda_xy + lambda_yx)
gamma_xzsp = Rz1p * (w_x + phi_x) + Rz2p * theta_x + Rz3p * lambda_x
gamma_yzsp = Rz1p * (w_y + phi_y) + Rz2p * theta_y + Rz3p * lambda_y
varepsilon_xxsp = 1 / 2 * gamma_xysp
varepsilon_xzsp = 1 / 2 * gamma_xzsp
varepsilon_yzsp = 1 / 2 * gamma_xzsp
########################## Strain on the bottom surface ###################################
varepsilon_xxsn = (Rn1 + h / 2) * w_xx + Rn1 * phi_xx + Rn2 * theta_xx + Rn3 * lambda_xx
varepsilon_yysn = (Rn1 + h / 2) * w_yy + Rn1 * phi_yy + Rn2 * theta_yy + Rn3 * lambda_yy
gamma_xysn = 2 * (Rn1 + h / 2) * w_xy + Rn1 * (phi_xy + phi_yx) + Rn2 * (theta_xy + theta_yx) + Rn3 * (lambda_xy + lambda_yx)
gamma_xzsn = Rz1n * (w_x + phi_x) + Rz2n * theta_x + Rz3n * lambda_x
gamma_yzsn = Rz1n * (w_y + phi_y) + Rz2n * theta_y + Rz3n * lambda_y
varepsilon_xysn = 1 / 2 * gamma_xysn
varepsilon_xzsn = 1 / 2 * gamma_xzsn
varepsilon_yzsn = 1 / 2 * gamma_yzsn
Mxx = integrate(sigma_xx * (R1 - z), (z, -h / 2, h / 2))
Myy = integrate(sigma_yy * (R1 - z), (z, -h / 2, h / 2))
Mxy = integrate(sigma_xy * (R1 - z), (z, -h / 2, h / 2))
Pxx1 = integrate(sigma_xx * R1, (z, -h / 2, h / 2))
Pxx2 = integrate(sigma_xx * R2, (z, -h / 2, h / 2))
Pxx3 = integrate(sigma_xx * R3, (z, -h / 2, h / 2))
Pyy1 = integrate(sigma_yy * R1, (z, -h / 2, h / 2))
Pyy2 = integrate(sigma_yy * R2, (z, -h / 2, h / 2))
Pyy3 = integrate(sigma_yy * R3, (z, -h / 2, h / 2))
Pxy1 = integrate(sigma_xy * R1, (z, -h / 2, h / 2))
Pxy2 = integrate(sigma_xy * R2, (z, -h / 2, h / 2))
Pxy3 = integrate(sigma_xy * R3, (z, -h / 2, h / 2))
Qx1 = integrate(kar * sigma_xz * Rz1, (z, -h / 2, h / 2))
Qx2 = integrate(kar * sigma_xz * Rz2, (z, -h / 2, h / 2))
Qx3 = integrate(kar * sigma_xz * Rz3, (z, -h / 2, h / 2))
Qy1 = integrate(kar * sigma_yz * Rz1, (z, -h / 2, h / 2))
Qy2 = integrate(kar * sigma_yz * Rz2, (z, -h / 2, h / 2))
Qy3 = integrate(kar * sigma_yz * Rz3, (z, -h / 2, h / 2))
Axx = Mxx.coeff(w_xx)
Bxx = Mxx.coeff(w_yy)
C1xx = Mxx.coeff(phi_xx)
C2xx = Mxx.coeff(theta_xx)
C3xx = Mxx.coeff(lambda_xx)
D1xx = Mxx.coeff(phi_yy)
D2xx = Mxx.coeff(theta_yy)
D3xx = Mxx.coeff(lambda_yy)
E1xx = Mxy.coeff(w_xy)
F1xx = Mxy.coeff(phi_xy)
F2xx = Mxy.coeff(theta_xy)
F3xx = Mxy.coeff(lambda_xy)
G1xx = Pxx1.coeff(w_xx)
H1xx = Pxx1.coeff(w_yy)
I11xx = Pxx1.coeff(phi_xx)
I12xx = Pxx1.coeff(theta_xx)
I13xx = Pxx1.coeff(lambda_xx)
J11xx = Pxx1.coeff(phi_yy)
J12xx = Pxx1.coeff(theta_yy)
J13xx = Pxx1.coeff(lambda_yy)
G2xx = Pxx2.coeff(w_xx)
H2xx = Pxx2.coeff(w_yy)
I21xx = Pxx2.coeff(phi_xx)
I22xx = Pxx2.coeff(theta_xx)
I23xx = Pxx2.coeff(lambda_xx)
J21xx = Pxx2.coeff(phi_yy)
J22xx = Pxx2.coeff(theta_yy)
J23xx = Pxx2.coeff(lambda_yy)
G3xx = Pxx3.coeff(w_xx)
H3xx = Pxx3.coeff(w_yy)
I31xx = Pxx3.coeff(phi_xx)
I32xx = Pxx3.coeff(theta_xx)
I33xx = Pxx3.coeff(lambda_xx)
J31xx = Pxx3.coeff(phi_yy)
J32xx = Pxx3.coeff(theta_yy)
J33xx = Pxx3.coeff(lambda_yy)
K1xx = Pxy1.coeff(w_xy)
L11xx = Pxy1.coeff(phi_xy)
L12xx = Pxy1.coeff(theta_xy)
L13xx = Pxy1.coeff(lambda_xy)
K2xx = Pxy2.coeff(w_xy)
L21xx = Pxy2.coeff(phi_xy)
L22xx = Pxy2.coeff(theta_xy)
L23xx = Pxy2.coeff(lambda_xy)
K3xx = Pxy3.coeff(w_xy)
L31xx = Pxy3.coeff(phi_xy)
L32xx = Pxy3.coeff(theta_xy)
L33xx = Pxy3.coeff(lambda_xy)
S1xx = Qx1.coeff(w_x)
S2xx = Qx2.coeff(w_x)
S3xx = Qx3.coeff(w_x)
T11xx = Qx1.coeff(phi_x)
T12xx = Qx1.coeff(theta_x)
T13xx = Qx1.coeff(lambda_x)
T21xx = Qx2.coeff(phi_x)
T22xx = Qx2.coeff(theta_x)
T23xx = Qx2.coeff(lambda_x)
T31xx = Qx3.coeff(phi_x)
T32xx = Qx3.coeff(theta_x)
T33xx = Qx3.coeff(lambda_x)
A11 = integrate2((I11xx * D3Xm * Yn + L11xx * D1Xm * D2Yn - T11xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A12 = integrate2((I12xx * D3Xm * Yn + L12xx * D1Xm * D2Yn - T12xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A13 = integrate2((I13xx * D3Xm * Yn + L13xx * D1Xm * D2Yn - T13xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A14 = integrate2(((J11xx + L11xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A15 = integrate2(((J12xx + L12xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A16 = integrate2(((J13xx + L13xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A17 = integrate2((G1xx * D3Xm * Yn + (H1xx + K1xx) * D1Xm * D2Yn - S1xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A21 = integrate2((I21xx * D3Xm * Yn + L21xx * D1Xm * D2Yn - T21xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A22 = integrate2((I22xx * D3Xm * Yn + L22xx * D1Xm * D2Yn - T22xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A23 = integrate2((I23xx * D3Xm * Yn + L23xx * D1Xm * D2Yn - T23xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A24 = integrate2(((J21xx + L21xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A25 = integrate2(((J22xx + L22xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A26 = integrate2(((J23xx + L23xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A27 = integrate2((G2xx * D3Xm * Yn + (H2xx + K2xx) * D1Xm * D2Yn - S2xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A31 = integrate2((I31xx * D3Xm * Yn + L31xx * D1Xm * D2Yn - T31xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A32 = integrate2((I32xx * D3Xm * Yn + L32xx * D1Xm * D2Yn - T32xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A33 = integrate2((I33xx * D3Xm * Yn + L33xx * D1Xm * D2Yn - T33xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A34 = integrate2(((J31xx + L31xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A35 = integrate2(((J32xx + L32xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A36 = integrate2(((J33xx + L33xx) * D1Xm * D2Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A37 = integrate2((G3xx * D3Xm * Yn + (H3xx + K3xx) * D1Xm * D2Yn - S3xx * D1Xm * Yn) * D1Xm * Yn, (x, 0, a), (y, 0, b))
A41 = integrate2(((J11xx + L11xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A42 = integrate2(((J12xx + L12xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A43 = integrate2(((J13xx + L13xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A44 = integrate2((I11xx * Xm * D3Yn + L11xx * D2Xm * D1Yn - T11xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A45 = integrate2((I12xx * Xm * D3Yn + L12xx * D2Xm * D1Yn - T12xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A46 = integrate2((I13xx * Xm * D3Yn + L13xx * D2Xm * D1Yn - T13xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A47 = integrate2((G1xx * Xm * D3Yn + (H1xx + K1xx) * D2Xm * D1Yn - S1xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A51 = integrate2(((J21xx + L21xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A52 = integrate2(((J22xx + L22xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A53 = integrate2(((J23xx + L23xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A54 = integrate2((I21xx * Xm * D3Yn + L21xx * D2Xm * D1Yn - T21xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A55 = integrate2((I22xx * Xm * D3Yn + L22xx * D2Xm * D1Yn - T22xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A56 = integrate2((I23xx * Xm * D3Yn + L23xx * D2Xm * D1Yn - T23xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A57 = integrate2((G2xx * Xm * D3Yn + (H2xx + K2xx) * D2Xm * D1Yn - S2xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A61 = integrate2(((J31xx + L31xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A62 = integrate2(((J32xx + L32xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A63 = integrate2(((J33xx + L33xx) * D2Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A64 = integrate2((I31xx * Xm * D3Yn + L31xx * D2Xm * D1Yn - T31xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A65 = integrate2((I32xx * Xm * D3Yn + L32xx * D2Xm * D1Yn - T32xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A66 = integrate2((I33xx * Xm * D3Yn + L33xx * D2Xm * D1Yn - T33xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A67 = integrate2((G3xx * Xm * D3Yn + (H3xx + K3xx) * D2Xm * D1Yn - S3xx * Xm * D1Yn) * Xm * D1Yn, (x, 0, a), (y, 0, b))
A71 = integrate2((C1xx * D4Xm * Yn + (D1xx + 2 * F1xx) * D2Xm * D2Yn - T11xx * D2Xm * Yn) * Xm * Yn, (x, 0, a), (y, 0, b))
A72 = integrate2((C2xx * D4Xm * Yn + (D2xx + 2 * F2xx) * D2Xm * D2Yn - T12xx * D2Xm * Yn) * Xm * Yn, (x, 0, a), (y, 0, b))
A73 = integrate2((C3xx * D4Xm * Yn + (D3xx + 2 * F3xx) * D2Xm * D2Yn - T13xx * D2Xm * Yn) * Xm * Yn, (x, 0, a), (y, 0, b))
A74 = integrate2((C1xx * Xm * D4Yn + (D1xx + 2 * F1xx) * D2Xm * D2Yn - T11xx * Xm * D2Yn) * Xm * Yn, (x, 0, a), (y, 0, b))
A75 = integrate2((C2xx * Xm * D4Yn + (D2xx + 2 * F2xx) * D2Xm * D2Yn - T12xx * Xm * D2Yn) * Xm * Yn, (x, 0, a), (y, 0, b))
A76 = integrate2((C3xx * Xm * D4Yn + (D3xx + 2 * F3xx) * D2Xm * D2Yn - T13xx * Xm * D2Yn) * Xm * Yn, (x, 0, a), (y, 0, b))
A77 = integrate2((Axx * (D4Xm * Yn + Xm * D4Yn) + 2 * (Bxx + E1xx) * D2Xm * D2Yn - S1xx * (Xm * D2Yn + D2Xm * Yn)) * Xm * Yn, (x, 0, a), (y, 0, b))
e1 = integrate2(D2Xm * Yn * Xm * Yn, (x, 0, a), (y, 0, b))
e2 = integrate2(Xm * D2Yn * Xm * Yn, (x, 0, a), (y, 0, b))
e3 = integrate2(D4Xm * Yn * Xm * Yn, (x, 0, a), (y, 0, b))
e4 = integrate2(D2Xm * D2Yn * Xm * Yn, (x, 0, a), (y, 0, b))
e5 = integrate2(Xm * D4Yn * Xm * Yn, (x, 0, a), (y, 0, b))
AP77 = (e1 + e2) * NCr
AK = np.array([[A77 + AP77]])
AK = Matrix(AK)
kk = solve(AK.det(), k)
print('Bulkling intensity factor for Kirchoff plate is k=', kk)
kv = kk[0]
print('The critical pressure is: Ncr=', (pi ** 2 * Diff / (a ** 2) * kv).evalf(6))参考文献^1: Tong L H, Lin F, Xiang Y, et al. Buckling analysis of nanoplates based on a generic third-order plate theory with shear-dependent non-isotropic surface stressesJ. Composite Structures, 2021, 265: 113708. https://doi.org/10.1016/j.compstruct.2021.113708
^2: 弹性力学：Rayleigh-Ritz法, 吃白饭的休伯利安号，https://www.eatrice.cn/post/RayleighRitzMethod/
选题思路理工科有着大量的数值计算的需求，现有的大部分的科学计算软件如matlab或mathmatica等均存在体积庞大、使用授权昂贵等问题。而python作为一款开源软件，其轻量、拓展性好、容易上手等完败那些难学的科学计算软件。同时python的用途广泛，学一门语言不仅可以做数值计算、还可以做数据挖掘、人工智能、其他工业软件插件开发等，对于非计算机科班出生的同学性价比极高。
本文介绍了python一款很受欢迎的符号计算模块：sympy，能够让读者了解python数值计算的优势，同时给出了常用功能的简单介绍，使得读者能够对python符号计算有一个完整且直观的理解。最后基于一篇论文的公式推导过程，给出了一个基于弹性力学的符号计算应用案例，更加直观地展现出python符号计算的强大以及其特别的魅力。
创作提纲为什么要使用python进行计算（分析当前常用方法的缺点，指出python计算的优点，引出sympy计算模块）sympy的安装与使用（介绍如何安装sympy）sympy的常用功能（通过高等数学和线性代数的常见计算场景介绍sympy，使得表达更加直观）sympy实际应用案例介绍（详细介绍了复杂公式的推导过程，并给出了相应的计算代码，展示将sympy投入实际应用的效果）参考文献（补充说明资料，数值计算往往是学科融合，需要一定的前置知识）