前往小程序，Get更优阅读体验！
立即前往
文档建议反馈控制台
首页
学习
活动
专区
工具
TVP
最新优惠活动
文章/答案/技术大牛
发布
首页
学习
活动
专区
工具
TVP最新优惠活动
返回腾讯云官网
社区首页 >专栏 >数值计算用Matlab？不，用python ｜ 技术创作特训营第一期

数值计算用Matlab？不，用python ｜ 技术创作特训营第一期

原创
作者头像
EatRice
发布于 2023-08-11 23:06:48
6800
发布于 2023-08-11 23:06:48
举报
文章被收录于专栏：个人学习笔记个人学习笔记

1 数值计算用什么

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

  1. Matlab是非开源的国外商业软件，存在安全问题以及盗版侵权问题；
  2. Matlab的安装包极大，对电脑的的要求很高；
  3. Matlab的跨平台能力较弱，编写出来的程序往往只能在安装了Matlab的机器上运行。
Matlab占用很高
Matlab占用很高

为了解决这些缺点，我们可以转而使用python来编写数值计算程序，当前的python版本支持多进程和多线程计算，numpy和sympy等高性能计算模块的开源共享使得python程序的计算性能和速度已经不输于matlab了。且python是开源的免费软件，一直为国内外的大神所维护能够保证性能和安全；python最新版的安装包才40M左右，十分轻量，在低配电脑上也有十分不错的兼容性；python可以在目前已知所有的操作系统上运行，编写一套代码可以在几乎所有的平台上运行。

基于以上的优点，这里强烈推荐一款python的计算模块：Sympy，他能够实现可视化的符号运算，并且与ipython兼容性十分不错，能够输出latex的可视化计算结果，如下图所示。本文将简要介绍Sympy的常用功能，并基于弹性力学给出一个计算模型作为算例，用于演示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 sympy

2.2 在jupyter notebook中显示公式

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

简单的latex公式
简单的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.359

3.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 polynomial

3.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
expr

3.8 符号化表达式

代码语言：python
复制
# 符号化表达式
str_expr='(x+1)**2'
expr=sp.sympify(str_expr)
expr

3.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_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$`表达式
`$\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符号计算的强大以及其特别的魅力。

创作提纲

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

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

2023腾讯·技术创作特训营 第一期

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

2023腾讯·技术创作特训营 第一期
#python
#sympy
#数值计算
#弹性力学
评论
登录后参与评论
0 条评论
热度
最新
登录 后参与评论
推荐阅读
LV.
文章
0
获赞
0
目录
  • 1 数值计算用什么
  • 2 sympy的安装与使用
    • 2.1 安装
      • 2.2 在jupyter notebook中显示公式
      • 3 sympy常用功能
        • 3.1 申明变量
          • 3.2 函数表达式（Expr）
            • 3.2.1 构造函数
            • 3.2.2 变量替换和数字赋值
            • 3.2.3 精确求值
            • 3.2.4微分
            • 3.2.5 积分
          • 3.3 极限
            • 3.4 解方程
              • 3.5 泰勒展开（不常见，但要会用）
                • 3.5.1 一元展开
                • 3.5.2 多元展开
                • 3.5.3 查看展开项的系数
              • 3.6 e的展开级数并化简
                • 3.6.1 e的指数级数的展开
                • 3.6.2 化简
              • 3.7 表达式具体输入值
                • 3.8 符号化表达式
                  • 3.9 极限
                    • 3.10 计算导数
                      • 3.10.1 对x求两次导
                      • 3.10.2 对多元函数求偏导
                    • 3.11 积分运算(integrate)
                      • 3.11.1 不定积分
                      • 3.11.2 定积分
                    • 3.12 解方程
                      • 3.13 解微分方程
                        • 3.14 矩阵运算
                          • 3.14.1 伴随矩阵
                      • 4 通过Rayleigh-Ritz法应用板理论计算板的变形
                      • 参考文献
                      • 选题思路
                      • 创作提纲
                      领券

                      腾讯云开发者

                      扫码关注腾讯云开发者

                      扫码关注腾讯云开发者

                      领取腾讯云代金券

                      热门产品

                      热门推荐

                      更多推荐

                      Copyright ? 2013 - 2024 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有?

                      深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059?深公网安备号 44030502008569

                      腾讯云计算（北京）有限责任公司 京ICP证150476号 | ?京ICP备11018762号 | 京公网安备号11010802020287

                      问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档

                      Copyright ? 2013 - 2024 Tencent Cloud.

                      All Rights Reserved. 腾讯云 版权所有

                      登录 后参与评论
                      http://www.vxiaotou.com