数值计算用Matlab?不,用python

1 数值计算用什么

作为理工科的社畜,懂计算会计算是一个必不可少的技能,其中尤其是对于土木工程人来说,结构力学、弹塑性力学、计算力学是数值计算中无法逾越的一道坎。由于Matlab简单使用,好学好操作,工科人往往都喜欢使用Matlab来实现数值算法。但是Matlab有几个缺点:

  1. Matlab是非开源的国外商业软件,存在安全问题以及盗版侵权问题;
  2. Matlab的安装包极大,对电脑的的要求很高;
  3. 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直接安装:

1
2
3
4
5
# 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文档:

1
2
3
4
5
6
7
8
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方法将字符串声明为符号变量,。

1
2
3
4
5
6
7
8
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构造抽象函数

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

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

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

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

3.2.4微分

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

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

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

3.3 极限

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

1
2
3
4
5
6
7
#####  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)

1
2
3
4
5
6
7
8
#### 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可以实现泰勒展开,具体函数抽象函数都可以。但是不能对多元函数同时泰勒展开。

1
2
3
4
5
6
7
8
#### 一元展开
# 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 多元展开

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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的展开级数并化简

1
2
3
4
5
6
7
8
# 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 符号化表达式

1
2
3
4
# 符号化表达式
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 计算导数

1
2
3
4
5
6
# 计算导数
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 解方程

1
2
3
4
5
6
7
8
9
10
11
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 解微分方程

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

1
2
3
4
5
6
7
8
9
10
11
12
#### 矩阵运算
# 构造矩阵
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 伴随矩阵

1
2
3
4
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$

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
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符号计算的强大以及其特别的魅力。

创作提纲

  1. 为什么要使用python进行计算(分析当前常用方法的缺点,指出python计算的优点,引出sympy计算模块)
  2. sympy的安装与使用(介绍如何安装sympy)
  3. sympy的常用功能(通过高等数学和线性代数的常见计算场景介绍sympy,使得表达更加直观)
  4. sympy实际应用案例介绍(详细介绍了复杂公式的推导过程,并给出了相应的计算代码,展示将sympy投入实际应用的效果)
  5. 与案例相关的参考文献(补充说明资料,数值计算往往是学科融合,需要一定的前置知识)
    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日
许可协议