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直接安装：

# Install the sympy modules using pip
pip install sympy

# Install the sympy modules using conda
conda install -c sympy

2.2 在jupyter notebook中显示公式

ipython的jupyter notebook支持加载mathjax脚本，能够实现可视化展示latex公式。在使用sympy可视化展示公式时，可以直接通过定义符号变量，并进行相关的运算来实现复杂公式的呈现，如下图所示：

简单的latex公式

当然也可以直接输出latex代码以嵌入至latex文档：

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方法将字符串声明为符号变量，。

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 函数表达式通过变量的运算构造具体函数，或者通过Function构造抽象函数

#### 函数表达式通过变量的运算构造具体函数，或者通过Function函数构造抽象函数。
# 具体函数
f=sympy.sqrt(3*x*y)+x*sympy.sin(y)+y**2+x**3
# 抽象函数
u=sympy.function('u')

3.2.2 expr.subs可以实现变量替换，替换成数字实现赋值

#### 变量替换与赋值
# 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位有效数字的精确值

#### 精确值
# expr.evalf(n)可以求一个表达式的保留n位有效数字的精确值
g3=f.subs({x:1,y:2})
print(g.evalf(4))  # 保留n位有效数字的精确值，8.359

3.2.4微分

sympy可以实现求微分，方法如下

### 微分
# 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可以实现自动求不定积分和定积分，区别在于是否传入积分上下限

#### 积分
# 可以实现自动求不定积分和定积分，区别在于是否传入积分上下限
l1=sympy.integrate(f,x)  #不定积分
l2=sympy.integrate(f,(x,1,3))    # 定积分

3.3 极限

sympy可以实现求极限，注意极限方向

#####  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)

#### 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可以实现泰勒展开，具体函数抽象函数都可以。但是不能对多元函数同时泰勒展开。

#### 一元展开
# 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 多元展开

函数的多元泰勒展开可以参考如下的代码。

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 polynomial

3.5.3 查看展开项的系数

1	`taylor.coeff(x) # 查看taylor1中项（x-x0)的系数`

3.6 e的展开级数并化简

# 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 表达式具体输入值

1
2
3

# 表达式输入具体值
expr=sp.exp(x)+1
expr

3.8 符号化表达式

# 符号化表达式
str_expr='(x+1)**2'
expr=sp.sympify(str_expr)
expr

3.9 极限

1
2
3

# 极限
sp.Sum(1/x**2,(x,1,sp.oo)).doit()############ 什么意思
sp.limit((1+1/x)**x,x,sp.oo)

由于结果显示最后一步，我拆开截屏（有的是因为在其上边拓展，有的是相同知识点的两个式子，没有将其分开，在编写代码运行的时候写在一个代码块中）

3.10 计算导数

# 计算导数
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 解方程

# 解方程
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 解微分方程

# 解微分方程
# 建立函数变量
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 矩阵运算

#### 矩阵运算
# 构造矩阵
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 伴随矩阵

A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
A
adj_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 $$

$$\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$

基于以上的推导过程，编写相应的计算程序求解：

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))

参考文献

选题思路

理工科有着大量的数值计算的需求，现有的大部分的科学计算软件如matlab或mathmatica等均存在体积庞大、使用授权昂贵等问题。而python作为一款开源软件，其轻量、拓展性好、容易上手等完败那些难学的科学计算软件。同时python的用途广泛，学一门语言不仅可以做数值计算、还可以做数据挖掘、人工智能、其他工业软件插件开发等，对于非计算机科班出生的同学性价比极高。

本文介绍了python一款很受欢迎的符号计算模块：sympy，能够让读者了解python数值计算的优势，同时给出了常用功能的简单介绍，使得读者能够对python符号计算有一个完整且直观的理解。最后基于一篇论文的公式推导过程，给出了一个基于弹性力学的符号计算应用案例，更加直观地展现出python符号计算的强大以及其特别的魅力。

创作提纲

为什么要使用python进行计算（分析当前常用方法的缺点，指出python计算的优点，引出sympy计算模块）
sympy的安装与使用（介绍如何安装sympy）
sympy的常用功能（通过高等数学和线性代数的常见计算场景介绍sympy，使得表达更加直观）
sympy实际应用案例介绍（详细介绍了复杂公式的推导过程，并给出了相应的计算代码，展示将sympy投入实际应用的效果）
与案例相关的参考文献（补充说明资料，数值计算往往是学科融合，需要一定的前置知识）
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 stresses[J]. Composite Structures, 2021, 265: 113708. https://doi.org/10.1016/j.compstruct.2021.113708 ↩
2. 弹性力学：Rayleigh-Ritz法, 吃白饭的休伯利安号，https://www.eatrice.cn/post/RayleighRitzMethod/ ↩

#学习 #内容分享 #数学

数值计算用Matlab？不，用python

https://www.eatrice.cn/post/MatlabOrPython/

作者

吃白饭舰长, 小张同学

发布于

2023年8月7日

许可协议

新的我和新的北京上一篇

弹性力学：Rayleigh-Ritz法下一篇