数值计算用Matlab？不，用python
1 数值计算用什么作为理工科的社畜，懂计算会计算是一个必不可少的技能，其中尤其是对于土木工程人来说，结构力学、弹塑性力学、计算力学是数值计算中无法逾越的一道坎。由于Matlab简单使用，好学好操作，工科人往往都喜欢使用Matlab来实现数值算法。但是Matlab有几个缺点：
Matlab是非开源的国外商业软件，存在安全问题以及盗版侵权问题；
Matlab的安装包极大，对电脑的的要求很高；
Matlab的跨平台能力较弱，编写出来的程序往往只能在安装了Matlab的机器上运行。
为了解决这些缺点，我们可以转而使用python来编写数值计算程序，当前的python版本支持多进程和多线程计算，numpy和sympy等高性能计算模块的开源共享使得python程序的计算性能和速度已经不输于matlab了。且python是开源的免费软件，一直为国内外的大神所维护能够保证性能和安全；python最新版的安装包才40M左右，十分轻量，在低配电脑上也有十分不错的兼容性；python可以在目前已知所有的操作系统上运行，编写一套代码可以在几乎所有的平台上运行。
基于以上的优点，这里强烈推荐一款python的计算模块：Sympy，他能够实现可视化的符号运算，并且与ipython兼容性十分不错，能够输出latex的可视化计算结果，如下图所示。本文将简要介绍Sympy的常用功能，并基于弹性力学给出一个计算模型作为算例，用于演示sympy在理工科的应用实战。
2 sympy的安装与使用sympy是一个开源模块，开源地址在github.com/sympy，代码包含详细的功能文档，建议直接fork下载查看。
2.1 安装sympy已经进入了PyPI，可以使用pip或conda直接安装：
12345# Install the sympy modules using pippip install sympy# Install the sympy modules using condaconda install -c sympy
2.2 在jupyter notebook中显示公式ipython的jupyter notebook支持加载mathjax脚本，能够实现可视化展示latex公式。在使用sympy可视化展示公式时，可以直接通过定义符号变量，并进行相关的运算来实现复杂公式的呈现，如下图所示：
当然也可以直接输出latex代码以嵌入至latex文档：
12345678from 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方法将字符串声明为符号变量，。
12345678import 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 函数表达式通过变量的运算构造具体函数，或者通过Function构造抽象函数12345#### 函数表达式通过变量的运算构造具体函数，或者通过Function函数构造抽象函数。# 具体函数f=sympy.sqrt(3*x*y)+x*sympy.sin(y)+y**2+x**3# 抽象函数u=sympy.function('u')
3.2.2 expr.subs可以实现变量替换，替换成数字实现赋值12345#### 变量替换与赋值# expr.subs()可以实现变量替换，替换成数字实现赋值。g1=f.subs(x,y) # 将f表达式中的x换成y,并将替换的结果赋给gg2=f.subs({x:2*x,y:2*y}) # 多次替换，字典g3=f.subs({x:1,y:2})
3.2.3 精确求值expr.evalf((n))可以求一个表达式的保留n位有效数字的精确值
1234#### 精确值# expr.evalf(n)可以求一个表达式的保留n位有效数字的精确值g3=f.subs({x:1,y:2})print(g.evalf(4)) # 保留n位有效数字的精确值，8.359
3.2.4微分sympy可以实现求微分，方法如下
123456### 微分# 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可以实现自动求不定积分和定积分，区别在于是否传入积分上下限
1234#### 积分# 可以实现自动求不定积分和定积分，区别在于是否传入积分上下限l1=sympy.integrate(f,x) #不定积分l2=sympy.integrate(f,(x,1,3)) # 定积分
3.3 极限sympy可以实现求极限，注意极限方向
1234567##### sympy可以实现求极限，注意极限方向# 趋于无穷lim1=sympy.limit(f,x,sympy.oo)# 趋于0，默认值dir="0",也就是趋于+0lim2=sympy.limit(f,x,0)# 趋于0，默认值dir="+"调整为dir="_",也就是趋于-0lim3=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)
12345678#### sympy可以实现解方程，方法是令Expr=0，所以在解方程时，要先构造宇哥#### 等于0的左端项。返回结果是一个列表，每一项是一个解，如果是方程组，解#### 解列表每一项是一个元组，元组对应位置是对应自变量的值func=f-3# 返回f=3时x的值sympy.solve(func,x) # x**2+y**2=1,x+y=1sympy.solve([x**2+y**2-1,x+y-1],[x,y])
3.5 泰勒展开（不常见，但要会用）3.5.1 一元展开sympy可以实现泰勒展开，具体函数抽象函数都可以。但是不能对多元函数同时泰勒展开。
12345678#### 一元展开# 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 多元展开函数的多元泰勒展开可以参考如下的代码。
12345678910111213141516171819202122232425262728def 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 polynomial
3.5.3 查看展开项的系数1taylor.coeff(x) # 查看taylor1中项（x-x0)的系数
3.6 e的展开级数并化简12345678# 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 表达式具体输入值123# 表达式输入具体值expr=sp.exp(x)+1expr
3.8 符号化表达式1234# 符号化表达式str_expr='(x+1)**2'expr=sp.sympify(str_expr)expr
3.9 极限123# 极限sp.Sum(1/x**2,(x,1,sp.oo)).doit()############ 什么意思sp.limit((1+1/x)**x,x,sp.oo)
由于结果显示最后一步，我拆开截屏（有的是因为在其上边拓展，有的是相同知识点的两个式子，没有将其分开，在编写代码运行的时候写在一个代码块中）
3.10 计算导数123456# 计算导数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 解方程123456789101112# 解方程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)E2sp.solveset(E2,x,domain=sp.Complexes)sp.solveset(E2,x,domain=sp.Reals)
3.13 解微分方程12345# 解微分方程# 建立函数变量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 矩阵运算123456789101112#### 矩阵运算# 构造矩阵sp.Matrix([[1,-1],[2,3],[3,4]])sp.Matrix([1,2,3])# 转置sp.Matrix([1,2,3]).TA=sp.Matrix([[1,2],[-2,1]])B=sp.Matrix([[3,4],[-1,2]]).TA+BA*BA**2B**2 # 得出结论：B转置后B**2结果也会转置
EA.tanspose() 为EA的转置矩阵EA.H 为EA的共轭转置矩阵
3.14.1 伴随矩阵1234A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])Aadj_A=A.adjugate()adj_A
4 通过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
$$
其中涉及到所有的变量计算表达式见参考文献[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$
基于以上的推导过程，编写相应的计算程序求解：
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349import numpy as npfrom sympy import *# Double integraldef integrate2(f, x, y): g = integrate(f, x) g = integrate(g, y) return g# Double differentialdef 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 platesST = 0# Define the condition of the boundaryBC = 3m = n = 1C11 = 263E9C12 = 154E9C44 = 127E9h = 0.8a = 10b = 10# h = 50E-9# b = a = 10 * hE = 25.5E9# E = (C11 - C12) * (C11 + 2 * C12) / (C11 + C12)v = 0.2# v = C12 / (C11 + C12)# Modified coefficient of the shear stressk = 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 = C44Diff = E * h ** 3 / 12 / (1 - v ** 2)DF = E * h ** 3NCr = k * pi ** 2 * Diff / (a ** 2)kv = 0print('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 ** 4if 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 ** 4PT = 1for i in range(0, 1, 2): kar = 1 if PT == 2: kar = 5 / 6R1 = R2 = R3 = 0Rz1 = Rz2 = Rz3 = 0Rp1 = Rp2 = Rp3 = 0Rn1 = Rn2 = Rn3 = 0Rz1p = Rz2p = Rz3p = 0Rz1n = Rz2n = Rz3n = 0print('Plate Theory: Kirchoff Plate')########################## strain ###################################varepsilon_xx = (R1 - z) * w_xx + R1 * phi_xx + R2 * theta_xx + R3 * lambda_xxvarepsilon_yy = (R1 - z) * w_yy + R1 * phi_yy + R2 * theta_yy + R3 * lambda_yygamma_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_xgamma_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_xysigma_xz = G * gamma_xzsigma_yz = G * gamma_yz########################## Strain on the top surface ###################################varepsilon_xxsp = (Rp1 - h / 2) * w_xx + Rp1 * phi_xx + Rp2 * theta_xx + Rp3 * lambda_xxvarepsilon_yysp = (Rp1 - h / 2) * w_yy + Rp1 * phi_yy + Rp2 * theta_yy + Rp3 * lambda_yygamma_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_xgamma_yzsp = Rz1p * (w_y + phi_y) + Rz2p * theta_y + Rz3p * lambda_yvarepsilon_xxsp = 1 / 2 * gamma_xyspvarepsilon_xzsp = 1 / 2 * gamma_xzspvarepsilon_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_xxvarepsilon_yysn = (Rn1 + h / 2) * w_yy + Rn1 * phi_yy + Rn2 * theta_yy + Rn3 * lambda_yygamma_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_xgamma_yzsn = Rz1n * (w_y + phi_y) + Rz2n * theta_y + Rz3n * lambda_yvarepsilon_xysn = 1 / 2 * gamma_xysnvarepsilon_xzsn = 1 / 2 * gamma_xzsnvarepsilon_yzsn = 1 / 2 * gamma_yzsnMxx = 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) * NCrAK = 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))
参考文献选题思路理工科有着大量的数值计算的需求，现有的大部分的科学计算软件如matlab或mathmatica等均存在体积庞大、使用授权昂贵等问题。而python作为一款开源软件，其轻量、拓展性好、容易上手等完败那些难学的科学计算软件。同时python的用途广泛，学一门语言不仅可以做数值计算、还可以做数据挖掘、人工智能、其他工业软件插件开发等，对于非计算机科班出生的同学性价比极高。
本文介绍了python一款很受欢迎的符号计算模块：sympy，能够让读者了解python数值计算的优势，同时给出了常用功能的简单介绍，使得读者能够对python符号计算有一个完整且直观的理解。最后基于一篇论文的公式推导过程，给出了一个基于弹性力学的符号计算应用案例，更加直观地展现出python符号计算的强大以及其特别的魅力。
创作提纲
为什么要使用python进行计算（分析当前常用方法的缺点，指出python计算的优点，引出sympy计算模块）
sympy的安装与使用（介绍如何安装sympy）
sympy的常用功能（通过高等数学和线性代数的常见计算场景介绍sympy，使得表达更加直观）
sympy实际应用案例介绍（详细介绍了复杂公式的推导过程，并给出了相应的计算代码，展示将sympy投入实际应用的效果）
与案例相关的参考文献（补充说明资料，数值计算往往是学科融合，需要一定的前置知识）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 stresses[J]. Composite Structures, 2021, 265: 113708. https://doi.org/10.1016/j.compstruct.2021.113708
↩
弹性力学：Rayleigh-Ritz法, 吃白饭的休伯利安号，https://www.eatrice.cn/post/RayleighRitzMethod/ ↩