Haopeng's Blog

Share equ happiness

三者存在的条件

  • 形心:目标具有几何外形;
  • 质心:目标具有几何外形、具有质量;
  • 重心:目标具有几何外形、具有质量、身处重力场。

质心

定义:质心是多质点系统的质量中心。 特点:若对该点施力,系统会沿着力的方向运动、不会旋转。 质心位置:质点的位置对质量取加权平均值,求得质心位置。

对于离散系统,三维物体的质心位置由下式求得:

\[ \bar{x} = \frac{\sum{x_k m_k}}{\sum{m_k}}, \bar{y} = \frac{\sum{y_k m_k}}{\sum{m_k}}, \bar{z} = \frac{\sum{z_k m_k}}{\sum{m_k}} \]

对于连续系统,上式推广为: \[ \bar{x} = \frac{\int_{\Omega}{x \rho} \rm{d}V}{\int_{\Omega}{ \rho} \rm{d}V}, \bar{y} = \frac{\int_{\Omega}{y \rho} \rm{d}V}{\int_{\Omega}{ \rho} \rm{d}V}, \bar{z} = \frac{\int_{\Omega}{z \rho} \rm{d}V}{\int_{\Omega}{ \rho} \rm{d}V} \]

其中,\(\rho\) 是物体的密度,为 \(x,y,z\) 的函数。上式为在物体所在空间 \(\Omega\) 进行体积分。

形心

定义:形心又被称为几何中心。形心是将该物体分成“矩”相等的两部分的所有超平面的交点。 形心位置:等价于取加权平均。

对于二维图形,形心位置由下式给出: \[ \bar{x} = \frac{\int{x f(x)}\rm{d}x}{\int{f(x)}\rm{d}x}, \bar{y} = \frac{\int{y f(y)}\rm{d}y}{\int{f(y)}\rm{d}y} \]

上式分子部分为 \(x\)\(y\) 的一阶矩。其中,\(f(x)\) 为几何图形的横坐标位于 \(x\) 点时在 \(y\) 方向上的长度;\(f(y)\) 为几何图形的纵坐标位于 \(y\) 点时在 \(x\) 方向上的长度。

上式的等价形式为: \[ \bar{x} = \frac{\int_D{x}\rm{d}S}{\int_D{ }\rm{d}S} = \frac{\int_D{x}\rm{d}S}{A} \]

上式积分表示在几何图形所在二维空间 \(D\) 进行面积分;\(A\) 表示二维几何图形的面积。(\(y\) 方向写法与之类似)

将该式推广至三维空间可得: \[ \bar{x} = \frac{\int_\Omega{x}\rm{d}V}{\int_\Omega{ }\rm{d}V} = \frac{\int_\Omega{x}\rm{d}V}{V} \]

对比发现,上式积分中如果加入密度项,则得到的质心的坐标。即,对于密度均匀的物体,其质心和形心重合

重心

定义:重心是重力作用的平均位置。 特点:各质点相对于重心的位置矢量乘上各质点的重力之和(合力矩)为零。

在均匀的重力场中,重心等同于质心。在非均匀的重力场中,质心和重心往往不在同一点。

数据类型

字符 str

  • 引号:在字符中显示单引号,则编程时用双引号将字符括起来;在字符中显示双引号,则编程时用单引号将字符括起来; 字符中的单引号:course = "He's finger" 字符串中的双引号:words = 'Python is for "beginners"'
  • 大小写course.upper()course.lower()course.title()
  • 查找替换course.find('xxx'),course.replace('xxx','aaa')
  • 字符长度:使用len()
  • f与花括号:可以被用来格式化输出字符(在字符中嵌入变量)。
    1
    2
    3
    first_name='John'
    last_name='Smith'
    msg = f"{first_name} [{last_name}] is a coder"

整数 int 浮点数 float 复数 complex

数据索引

  1. Python中的索引从0开始
  2. 数据索引中的负号表示从结束端开始:[1:-1]是正确的,-1表示倒数第二个数据

数据切片索引 与Matlab不同的是,A[1:5]返回的不是数列A中位于1、2、3、4、5号位置的数据;取而代之的是:

1
2
3
4
5
6
7
8
9
10
A = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
print(A[1:5])

Output:
[2, 3, 4, 5]

-------------------------------------------
[A]: | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
cut: -0---1---2---3---4---5---6---7---8---9----10-

if

常用调用模式如下,用缩进控制开始和结束(与matlab不同):

1
2
3
4
5
6
if is_hot:
print('xxxxx')
elif is_coldL:
print('xxxxxx')
else:
print('xxxx')

常用逻辑运算符 1. 和:and,例如:is_hot and is_cold 1. 或:or 1. 非:not,例如:is_cold and not is_windnot会将其后的False转化为True;上例中如果is_cold = Trueis_wind = False,则返回True

比大小: 1. 大于:>、小于:<、等于:== 1. 大于等于:>=、小于等于:<= 1. 不等于:!=

循环

while常用调用模式如下,break可以用来跳出程序:

1
2
3
4
5
i=1
while i<=5:
print('i')
i = i + 1
print('Done') #用缩进控制while的结束

for常用调用模式

1
2
for item in [1,2,3,4]
print(item)

函数

常用函数

绝对值:abs() 四舍五入到整数:round() 生成等间距数列:range(5,10,2);step=2,from 5 to 10,输出 5 7 9

常用命令

运算符号

  1. +;减-;乘*;除/;幂**
  2. 特殊除号:得到结果的整数部分//;得到结果的余数%
  3. 增广运算符X += 3x = x + 3等价、x -= 3同理

import

调用包:import XXX,例如import math;使用时math.acos()

Questions

  1. methods 和 function 在Python有何区别的?

概述

加权残值伽辽金法 (Galerkin method of weighted residuals) 也被称作加权余量伽辽金法、伽辽金法 (Galerkin method),常被看作是加权残值(余量)法的一种。该方法是俄罗斯数学家鲍里斯·格里戈里耶维奇·伽辽金 (Boris Galerkin) 发明的一种数值分析方法。应用这种方法可以将求解微分方程问题简化成为线性方程组的求解问题。

伽辽金法原理概述: 采用微分方程对应的弱形式,通过选取有限多项试函数(又称基函数或形函数),将它们叠加;再要求结果在求解域内及边界上的加权积分(权函数为试函数本身)满足原方程,便可以将求解微分方程近似解的问题转化为求解一组线性代数方程,且自然边界条件能够自动满足。

必须强调指出的是,伽辽金法所得到的只是在原求解域内的一个近似解(仅仅是加权平均满足原方程,并非在每个点上都满足)。

为了阐明伽辽金法的原理及特点,本文先介绍相关概念——微分方程的等效积分形式弱形式,再介绍加权残值伽辽金法的基本步骤和特点。

等效积分弱形式

数学中, 微分方程的弱解或广义解是指对该方程中的微分可能不存在, 但是在某种精确定义的意义下满足该方程的解。 对于不同种类的微分方程, 弱解的定义性质也可能不同。

设函数 \(\mathbf{u}\) 在域 \(\Omega\) 中满足如下微分方程组:

\[ \mathbf{A}(\mathbf{u}) = 0 \]

\(\Omega\) 可能是体积域或面积域。同时未知函数 \(\mathbf{u}\) 还要在边界 \(\Gamma\) 满足边界条件:

\[ \mathbf{B}(\mathbf{u}) = 0 \]

由于微分方程组在域 \(\Omega\) 中的每一点都必须为0,因此有:

\[ \int _\Omega \mathbf{v} ^\mathrm{T} \mathbf{A}(\mathbf{u}) \mathrm{d} \Omega \equiv \int _\Omega \left( v _1 A _1 (\mathbf{u})+v _2 A _2 (\mathbf{u})+\cdots \right)\mathrm{d} \Omega \equiv 0 \]

上式的数字下标代表对应向量的分量,\(\mathbf{v}\) 是一组与微分方程个数相等的任意函数。

假如边界条件在边界上每一点也同时得到满足,则对任意一组函数 \(\bar{\mathbf{v}}\),下式成立:

\[ \int _\Gamma \bar{\mathbf{v}} ^\mathrm{T} \mathbf{B}(\mathbf{u}) \mathrm{d} \Gamma \equiv \int _\Gamma \left( \bar{v} _1 B _1 (\mathbf{u})+ \bar{v} _2 B _2 (\mathbf{u})+\cdots \right)\mathrm{d} \Gamma \equiv 0 \]

因此有:

\[ \int _\Omega \mathbf{v} ^\mathrm{T} \mathbf{A}(\mathbf{u}) \mathrm{d} \Omega + \int _\Gamma \bar{\mathbf{v}} ^\mathrm{T} \mathbf{B}(\mathbf{u}) \mathrm{d} \Gamma \equiv 0 \]

上式是与微分方程组 \(\mathbf{A}(\mathbf{u})=0\) 完全等效的等效积分形式

在很多情况下可以对上式进行分部积分得到另一种形式:

\[ \int _\Omega \mathbf{C} ^\mathrm{T} (\mathbf{v}) \mathbf{D}(\mathbf{u}) \mathrm{d} \Omega + \int _\Gamma \mathbf{E} ^\mathrm{T} (\bar{\mathbf{v}}) \mathbf{F}(\mathbf{u}) \mathrm{d} \Gamma \equiv 0 \]

其中 \(\mathbf{C}, \mathbf{D}, \mathbf{E}, \mathbf{F}\) 是微分算子,与分部积分之前的 \(\mathbf{A}\) 相比,他们所包含的导数的阶数较低,这样对函数 \(\mathbf{u}\) 只需求较低阶的连续性就可以了。上式中降低 \(\mathbf{u}\) 的连续性是以提高 \(\mathbf{v}\)\(\mathbf{v^{'}}\) 的连续性为代价的。然而,适当提高 \(\mathbf{v}\) 的连续性并不困难,因为它们是可以选择的已知函数。

这种通过适当提高对任意函数 \(\mathbf{v}\)\(\mathbf{v^{'}}\) 的连续性要求,以降低对微分方程场函数 \(\mathbf{u}\) 的连续性要求所建立的等效积分形式称为微分方程 \(\mathbf{A}\)等效积分弱形式。等效积分弱形式的解称为微分方程的弱解

加权残值伽辽金法

概述

采用使余量的加权积分为零来求得微分方程近似解的方法称为加权残值法,也称加权余量法。加权残值法是一种基于等效积分形式的近似方法。根据所选取权函数的不同,加权残值法可分为:配点法、子域法、最小二乘法、力矩法和加权残值伽辽金法

在应用数学中,加权残值法被用来求解微分方程的近似解,具体步骤如下: 1. 假设近似解的形式是有限个试探函数(或基函数、形函数) \(\mathbf{N} _i\) 的叠加,叠加时的系数为 \(\mathbf{a} _i\)。 1. 由于近似解不能精确满足原微分方程和边界条件,近似解带入原微分方程和边界条件后将产生残差/余量(等式左端不为零的项)。 1. 将得到的关于残差的方程写作等效积分形式,并用规定的权函数 \(\mathbf{W} _j\) 代替任意函数 \(\bar{v}\)。 1. 令得到的由残差和权函数表示的微分方程等于零,可求得待定系数 \(\mathbf{a} _i\),最终得到了原微分方程的近似解。

Notation:用近似解的试探函数 \(\mathbf{N} _j\) 作为权函数 \(\mathbf{W} _j\) 时 (\(\mathbf{W} _j = \mathbf{N} _j\)),上述步骤被称为加权残值伽辽金法。

具体推导

在求解域 \(\Omega\) 中,若场函数 \(\mathbf{u}\) 是精确解,则在域 \(\Omega\) 中的任何一点都满足微分方程 \(\mathbf{A}(\mathbf{u})=0\),同时在边界 \(\Gamma\) 上满足边界条件 \(\mathbf{B}(\mathbf{u})=0\),此时等效积分式或其弱形式必然也得到严格满足。但是对于复杂的实际问题,这样的精确解往往很难找到,因此需要设法找到具有一定精度的近似解。

假设未知场函数 \(\mathbf{u}\) 可以用近似函数表示。近似函数是带有待定参数的已知函数一般形式是:

\[ \mathbf{u} \approx \mathbf{u}^* = \sum ^n _{i=1} \mathbf{N} _i \mathbf{a} _i = \mathbf{N} \mathbf{a} \]

其中 \(\mathbf{a} _i\) 是待定参数,\(\mathbf{N} _i\) 是被称为试探函数(或基函数、形函数)的已知函数,它取自完全的函数序列,是线性独立的。近似解通常选择使之满足强制边界条件和连续性的要求。

完全函数序列:任一函数都可以用此序列表示。

显然,试探函数为有限项时,近似解是不能精确满足微分方程和边界条件的,它们将产生残差 \(\mathbf{R}\)\(\bar{\mathbf{R}}\)

\[ \mathbf{A}(\mathbf{Na}) = \mathbf{R} \]

\[ \mathbf{B}(\mathbf{Na}) = \bar{\mathbf{R}} \]

残差 \(\mathbf{R}\)\(\bar{\mathbf{R}}\) 也被称为余量。将原方程写作等效积分形式,并用 \(n\) 个规定的函数代替任意函数 \(\mathbf{v}\)\(\bar{\mathbf{v}}\)

\[ \mathbf{v}=\mathbf{W} _j, \quad \bar{\mathbf{v}}=\bar{\mathbf{W}} _j \quad (j=1,2,\dots,n) \]

则原微分方程的等效形式变为:

\[ \int _\Omega \mathbf{W} _j ^\mathrm{T} \mathbf{A}(\mathbf{Na}) \mathrm{d} \Omega + \int _\Gamma \bar{\mathbf{W}} _j ^\mathrm{T} \mathbf{B}(\mathbf{Na}) \mathrm{d} \Gamma = 0 \]

亦可写为余量形式:

\[ \int _\Omega \mathbf{W} _j ^\mathrm{T} \mathbf{R} \mathrm{d} \Omega + \int _\Gamma \bar{\mathbf{W}} _j ^\mathrm{T} \bar{\mathbf{R}} \mathrm{d} \Gamma = 0 \]

也可写出等效积分弱形式:

\[ \int _\Omega \mathbf{C} ^\mathrm{T} (\mathbf{W} _j) \mathbf{D}(\mathbf{Na}) \mathrm{d} \Omega + \int _\Gamma \mathbf{E} ^\mathrm{T} (\bar{\mathbf{W}} _j) \mathbf{F}(\mathbf{Na}) \mathrm{d} \Gamma \equiv 0 \]

以上三式的意义是:通过选择待定系数 \(\mathbf{a} _i\) ,强迫余量在某种平均意义上等于零。\(\mathbf{W} _j\)\(\bar{\mathbf{W}} _j\) 称为权函数。 余量的加权积分为零就得到了一组求解方程,用以求解近似解的待定系数 \(\mathbf{a}\),从而得到原问题的近似解答。

近似函数所取的试探函数的项数越多,近似解的精度越高。当项数趋于无穷时,近似解将收敛于精确解。

\(\mathbf{W} _j= \mathbf{N} _j\),在边界上 \(\bar{\mathbf{W}} _j= -\mathbf{W} _j = -\mathbf{N} _j\)。即简单的利用近似解的试探函数序列作为权函数时,以上步骤称为加权残值伽辽金法。此时有:

\[ \int _\Omega \mathbf{N} _j ^\mathrm{T} \mathbf{A}\left( \sum ^n _{i=1}\mathbf{N} _i \mathbf{a} _i \right) \mathrm{d} \Omega + \int _\Gamma \bar{\mathbf{N}} _j ^\mathrm{T} \mathbf{B}\left( \sum ^n _{i=1}\mathbf{N} _i \mathbf{a} _i \right) \mathrm{d} \Gamma = 0 \]

加权残值伽辽金法的特点

特点:如果算子 \(\mathbf{A}\) 是线性自伴随的,则采用伽辽金法得到的求解方程的系数矩阵是对称的。 在高维问题中,这种对称性会大大减少计算量。因此使用加权残值法建立有限元格式时,几乎毫无例外地采用伽辽金法。现以一维热传导问题为例说明该特点。

在一维热传导问题中,如果热传导系数取1,则微分方程为:

\[ A(\phi) = \frac{\mathrm{d}^2 \phi}{\mathrm{d} x^2}+Q(x) = 0 \quad (0 \leq x \leq L) \]

其中,当 \(0\leq x\leq \frac{L}{2}\)\(Q(x)=1\) ;当 \(\frac{L}{2}\leq x\leq L\)\(Q(x)=0\)。 边界条件为 \(x=0\)\(x=L\) 时,\(\phi =0\)

取傅里叶级数作为近似解:

\[ \phi \approx \tilde{\phi} = \sum ^n _{i=1} a _i \sin \frac{i\pi x}{L} \]

其中 \(a_i\) 为待定参数,试探函数 \(N_i= \sin \frac{i\pi x}{L}\)。近似解满足边界条件,因此在边界上不产生余量。将近似解带入原方程的等效积分形式中:

\[ \int ^L _0 W _j \left[ \frac{\mathrm{d}^2}{\mathrm{d} x^2} \left( \sum ^n _{i=1} N _i a _i \right) +Q \right] \mathrm{d} x=0 \quad (j=1,2,\dots,n) \]

采用伽辽金法时,因为权函数 \(W _j = N _j\) 是连续的,并且在两端等于 \(N _j =0\)。对上式分部积分可以得到等效积分形式的弱形式:

\[ \int ^L _0 \left[ \frac{\mathrm{d}W _j}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d} x} \left( \sum ^n _{i=1} N _i a _i \right) + W _j Q \right] \mathrm{d} x=0 \quad (j=1,2,\dots,n) \]

上式可以改写为:

\[ \mathbf{Ka-P=0} \]

其中,\(\mathbf{P}\)\(\mathbf{a}\)\(n\) 维向量

\[ P _j =\int ^L _0 W _j Q \mathrm{d}x \]

\[ K _{ij} = \frac{\mathrm{d}W _i}{\mathrm{d}x} \frac{\mathrm{d} N _j}{\mathrm{d} x} \]

可以看到,当 \(W _i= N _i\) 时,将有 \(K _{ij} = K _{ji}\),换而言之,采用伽辽金法得到的求解待定参数 \(a _i\) 的代数方程组的系数矩阵 \(\mathbf{K}\) 是对称的。

概述

最小势能原理和最小余能原理是物理和工程学中的基本概念,被应用于力学、电学、生物学等多种学科。

弹性力学有两种提法,一种微分提法,一种变分提法。微分提法从微元入手,得到严格满足方程和边界条件的强解;变分提法从整体(能量)入手,得到弱解。最小势能原理和最小余能原理是为后者服务的,它们从能量入手,为一些经典方法的提出奠定了理论基础,如Galerkin(伽辽金)方法、Rayleigh–Ritz(瑞利-里兹)方法等。本文只讨论弹性力学中的最小势能原理与最小余能原理。

为了阐明这两种原理及其关系,需要引入一些概念:运动许可状态、静力许可状态、虚位移、虚力、系统的总势能、系统的总余能、应变能密度、余应变能密度、可能功原理

相关概念的引入

基本概念

  • 运动许可状态: 满足位移约束条件的(可能存在的)状态。相应的变分是虚位移。
  • 静力许可状态: 满足平衡方程和力边界条件的(可能存在的)状态。相应的变分是虚力。
  • 虚位移: 给定的瞬时和位形上,虚位移是符合约束条件的无穷小位移。
  • 虚力: 给定的瞬时和位形上,虚力是符合平衡方程和力边界条件的的无穷小力。
  • 总势能 \(\Pi\) 等于应变能+外力势(\(\Pi=U+W\)
  • 总余能 \(\Pi^*\) 等于余应变能+外力余势(\(\Pi ^ * = U ^ *+ W ^ *\)

应变能、余应变能、外力势、外力余势的概念见下一节

应变能和余应变能

应变能: 弹性体受外力作用产生应变,产生应力场和位移场。外力对弹性体所做的功以变形能的形式储存在弹性体内,这种变形能通常称为应变能。

现以杆件的拉伸问题为例介绍应变能余应变能的概念。如下图(a)所示,杆件原长为 \(l\),伸长量为 \(u=\varepsilon l\),受到一个平行于杆件的静力 \(P=\sigma f\)\(f\) 为横截面积。

此时,外力 \(P\) 所做的功所引起的势能损失为图(b)中曲线下半部分面积 \(W\),被称为外力势(外载荷的势能):

\[ W= -\int _0^u P \mathrm{d} u = -fl \int _0^\varepsilon \sigma \mathrm{d} \varepsilon \]

上式中的正负号表示外力所做的功引起的系统势能的增减,此时外力方向与位移方向相同,外力做正功,系统有势能损失,为负号。这可能不好理解,那么我们考虑重力对物体的影响。当一个重物被抬升一段距离,此时系统势能增加,重力做负功。即外力方向与位移方向相反,外力做负功,系统势能增加。

画出系统的应力应变曲线图(c),则上式左端

\[ A = \int _0^\varepsilon \sigma \mathrm{d} \varepsilon \]

代表单位体积下的应变势能,定义为应变能密度。则系统的应变能 \(U\) 为:

\[ U= \int _V A \mathrm{d}V = fl \int _0^\varepsilon \sigma \mathrm{d} \varepsilon \]

与应变能相对的是余应变能

应变能 \(U\) 余应变能 \(U^*\)
应变能密度 \(A\) 余应变能密度 \(B\)
外力势 \(W\) 外力余势 \(W^*\)

余应变能密度 \(B\) 对应图(c)上半部分面积:

\[ B = \int _0^\sigma \varepsilon \mathrm{d} \sigma \]

则系统的余应变能 \(U^*\) 为:

\[ U= \int _V B \mathrm{d}V = fl \int _0^\sigma \varepsilon \mathrm{d} \sigma \]

外力余势 \(W^*\) 对应图(b)上半部分面积:

\[ W^* = -\int _0^P u \mathrm{d} P \]

外力余势通常在支撑部分。

应变能与余应变能之间的关系:

\[ A+B = \int _V \sigma \varepsilon \mathrm{d} V \]

可能功原理

考虑弹性体的两个状态,状态一是应力许可的,状态二是运动许可的。

利用上标 \(s\) 表示状态一(\(f _i\) 为体力,\(p _i\) 为面力)

\[ \sigma _{ij,j}^{(s)}+ f _i^{(s)}=0 \]

\[ \sigma _{ij}^{(s)} v _j =p _i^{(s)} \]

利用上标 \(k\) 表示状态二

\[ \varepsilon _{ij}^{(k)} = \frac{1}{2} \left( u _{i,j}^{(k)} + u _{j,i}^{(k)} \right) \]

考虑状态一中的体力、面力和可能应力在状态二的相应可能位移和可能应变上所做的功,分别为:

\[ A _f = \int _V f _i^{(s)} u _i^{(k)} \mathrm{d} V \]

\[ A _p = \int _S p _i^{(s)} u _i^{(k)} \mathrm{d} S = \int _S \sigma _{ij}^{(s)} u _i^{(k)} v _j \mathrm{d} S = \int _V \left( \sigma _{ij}^{(s)} u _i^{(k)} \right) _{,j}\mathrm{d} V = \int _V \sigma _{ij,j}^{(s)} u _i^{(k)} \mathrm{d} V +\int _V \sigma _{ij}^{(s)} u _{i,j}^{(k)} \mathrm{d} V \]

\[ A _\sigma = \int _V \sigma _{ij}^{(s)} \varepsilon _{ij}^{(k)} \mathrm{d} V \]

可能功原理: 外力功等于内力功,状态一的外力在状态二位移做的功等于状态一的应力在状态二的应变上做的功。

物理意义: 外力与内力构成自平衡力系,他们在可能位移上做的总功为零。 适用性: 静力可能状态和几何可能状态时完全独立的;与本构关系无关,适用于任何连续介质

推导过程如下:

\[ A _p = \int _V \sigma _{ij,j}^{(s)} u _i^{(k)} \mathrm{d} V +\int _V \sigma _{ij}^{(s)} u _{i,j}^{(k)} \mathrm{d} V = \int _V f _i^{(s)} u _i^{(k)} \mathrm{d} V + \int _V \sigma _{ij}^{(s)} \varepsilon _{ij}^{(k)} \mathrm{d} V = -A _f +A _\sigma \]

\[ A _f + A _p = A _\sigma \]

\[ \int _V f _i^{(s)} u _i^{(k)} \mathrm{d} V + \int _S p _i^{(s)} u _i^{(k)} \mathrm {d} S = \int _V \sigma _{ij}^{(s)} \varepsilon _{ij}^{(k)} \mathrm{d} V \]

最小势能原理

原理简述

最小势能原理: 在满足弹性体的几何方程和位移边界条件的所有容许的位移中,真实的位移必使弹性体的总势能有驻值。当弹性体为稳定平衡时,其总势能为极小值。

简述为:在一切运动许可状态中,真实状态的总势能最小。

原理推导

图中是一个一般意义下的弹性体。受到的外力分为面力 \(\bar{p} _i\) 和体力 \(f _i\);弹性体的外表面 \(S _u\) 所在区域受到位形约束 \(\bar{u} _i\)

\[ \varepsilon _{i j}=\frac{1}{2}\left(u _{i, j}+u _{j, i}\right) \]

\[ \sigma _{i j, j}+f _{i}=0 \]

\[ \sigma _{i j} v _{i}=\bar{p} _{j} \quad \text { on } S _{\sigma} \]

\[ u _{i}= \bar{u} _{i} \quad \text { on } S _{u} \]

\[ \sigma _{i j}=\frac{\partial A\left(\varepsilon _{i j}\right)}{\partial \varepsilon _{i j}} \]

真实状态的总势能 \(\Pi\) 分为应变能 \(U\) 和 外力势 \(W\) 两部分:

\[ \Pi = U+W = \int _{V} A \left( \varepsilon _{i j} \right) \mathrm{d} V- \left( \int _{V} f _{i} u _{i} \mathrm{d} V+\int _{S _{\sigma}} \bar{p} _{i} u _{i} \mathrm{d} S \right) \]

其中 \(A(\varepsilon_{ij})\) 为应变能密度。

运动许可状态\(\varepsilon _{ij} ^{(k)}\)\(u _{i} ^{(k)}\) 描述,仅满足变形关系:

\[ \varepsilon_{i j}^{(k)}=\frac{1}{2}\left(u_{i, j}^{(k)}+u_{j, t}^{(k)}\right), \quad u_{i}^{(k)}=\bar{u}_{i} \]

运动许可状态的总势能为:

\[ \Pi^{(k)}=\int _{V} A^{(k)}\left(\varepsilon _{i j}^{(k)}\right) \mathrm{d} V-\int _{V} f _{i} u _{i}^{(k)} \mathrm{d} V-\int _{S _{\sigma}} \bar{p} _{i} u _{i}^{(k)} \mathrm{d} S \]

两种状态的总势能差为:

\[ \begin{aligned} \Pi^{(k)}-\Pi &=\int_{V}\left[A^{(k)}\left(\varepsilon_{i j}^{(k)}\right)-A\left(\varepsilon_{i j}\right)\right] \mathrm{d} V \\ &-\int_{V} f_{i}\left[u_{i}^{(k)}-u_{i}\right] \mathrm{d} V-\int_{S} p_{i}\left[u_{i}^{(k)}-u_{i}\right] \mathrm{d} S \end{aligned} \]

注意到与面力相关项的积分区域 \(S_\sigma \rightarrow S\),这是因为两种状态都满足位形约束,所以两种状态在 \(S _u\) 的位移差恒为0。

利用可能功原理,外力项可以被改写为:

\[ -\int_{V} f_{i}\left[u_{i}^{(k)}-u_{i}\right] \mathrm{d} V-\int_{S} p_{i}\left[u_{i}^{(k)}-u_{i}\right] \mathrm{d} S =-\int_{V} \sigma_{i j}\left(\varepsilon_{i j}^{(k)}-\varepsilon_{i j}\right) \mathrm{d} V=-\int_{V} \frac{\partial A}{\partial \varepsilon_{i j}}\left(\varepsilon_{i j}^{(k)}-\varepsilon_{i j}\right) \mathrm{d} V \]

则两种状态的势能差为:

\[ \Pi^{(k)}-\Pi=\int_{V}\left[A^{(k)}\left(\varepsilon_{i j}^{(k)}\right)-A\left(\varepsilon_{i j}\right) - \frac{\partial A}{\partial \varepsilon_{i j}}\left(\varepsilon_{i j}^{(k)}-\varepsilon_{i j}\right)\right] \mathrm{d} V \]

只要应变能密度函数 \(A\) 是凸函数,则必有:

\[ \Pi^{(k)}-\Pi \geq 0 \]

最小势能原理的变分形式:

\[ \delta \Pi=0 \]

\[ \delta \Pi=\delta U+\delta V=\delta \int _{V} A \mathrm{d} V-\int _{V} f _{i} \delta u _{i} \mathrm{d} V-\int _{S _{\sigma}} \bar{p} _{i} \delta u _{i} \mathrm{d} S \]

\[ \delta U=-\delta V=\int _{V} f _{i} \delta u _{i} \mathrm{d} V+\int _{S _{\sigma}} \bar{p} _{i} \delta u _{i} \mathrm{d} S \]

最小余能原理

原理简述

最小势能原理: 在弹性体内满足平衡微分方程,在应力边界上满足应力边界条件的所有容许的应力状态中,真实的应力状态(满足变形协调条件的应力)必使其总余能为极小值。

简述为:在一切静力许可状态中,真实状态的总余能最小。

原理推导

还是考虑上图中的弹性体,应变可以用余应变能 \(B\) 表示:

\[ \varepsilon_{i j}=\frac{\partial B\left(\sigma_{i j}\right)}{\partial \sigma_{i j}} \]

真实状态的系统总余能 \(\Pi^*\) 由 余应变能 \(U^*\) 和外力余势 \(W^*\) 构成:

\[ \Pi ^*= U^* +W^* \int_{V} B \left(\sigma_{i j}\right) \mathrm{d} V- \int_{S_{u}} p_{i} \bar{u}_{i} \mathrm{~d} S \]

静力许可状态仅满足静力关系:

\[ \sigma _{i j, j}^{(s)}+f _{i}=0 \quad \text{ on } V \]

\[ p_{i}^{(s)}=\sigma_{i j}^{(s)} v_{j} \quad \text{ on } S _u \]

\[ p_{i}^{(s)}=\bar{p}_{i}\quad \text{ on } S _\sigma \]

静力许可状态的总余能为:

\[ \Pi ^{*(s)}= \int_{V} B^{(s)}\left(\sigma_{i j}^{(s)}\right) \mathrm{d} V- \int_{S_{u}} p_{i}^{(s)} \bar{u}_{i} \mathrm{~d} S \]

两种状态的势能差为:

\[ \Pi^{* (s)}-\Pi^*= \int _{V}\left[B^{(s)} \left( \sigma _{ij}^{(s)} \right)-B\left(\sigma _{i j}\right)\right] \mathrm{d} V -\int _{V}\left(f _{i}^{(s)}-f _{i}\right) u _{i} \mathrm{d} V-\int _{S}\left(p _{i}^{(s)}-p _{i}\right) u _{i} \mathrm{d} S \]

利用可能功原理,外力项可以被改写为:

\[ -\int _{V} \left( f _{i}^{(s)}-f _{i}\right) u _{i} \mathrm {d} V-\int _{S}\left(p _{i}^{(s)}-p _{i}\right) u _{i} \mathrm {d} S=-\int _{V} \varepsilon _{i j} \left( \sigma _{i j}^{(s)}-\sigma _{i j}\right) \mathrm {d} V =-\int _{V} \frac {\partial B}{\partial \sigma _{i j}}\left(\sigma _{i j}^{(s)}-\sigma _{i j}\right) \mathrm {d} V \]

只要应变余能密度函数 \(B\) 是凸函数,则必有:

\[ \Pi^{* (s)}-\Pi^* \geq 0 \]

最小势能原理和最小余能原理的区别与联系

联系:

\[ \Pi+\Pi^* =\int _{V} \left[ A+B \right] \mathrm{d} V-\int _{V} f _{i} u _{i} \mathrm{d} V-\int _{S _{\sigma}} \bar{p} _{i} u _{i} \mathrm{d} S-\int _{S _{u}} p _{i} \bar{u} _{i} \mathrm{d} S =\int _{V} \sigma _{i j} \varepsilon _{i j} \mathrm{d} V-\int _{V} f _{i} u _{i} \mathrm{d} V-\int _{S} p _{i} u _{i} \mathrm{d} S \]

由可能功原理可得 \(\Pi +\Pi^*=0\)

于是有:

\[ \Pi^{* (s)} \geq \Pi^* =-\Pi \geq - \Pi^{(k)} \]

因此用这两种原理建立数值方程求解有如下特点: * 利用最小势能原理求得位移近似解的弹性变形能是精确解变形能的下届,即近似的位移场在总体上偏小,也就是说结构的计算模型显得偏于刚硬。 * 利用最小余能原理得到的近似解的弹性余能是精确余能的上界,即近似的应力解在总体上偏大,结构计算模型模型偏于柔软。 * 分别利用这两个极值原理求解同一个问题时,可以获得这个问题的上界和下界,可以较准确的估算所得近似解的误差。

区别: 最小势能原理以位移为基本变量,要求位移场事先满足几何方程和给定位移边界条件。最小余能原理以应力为基本变量,要求应力场事先满足平衡方程和给定面力边界条件。

用这两个原理求解的优点是通常只有一个场函数,且泛函具有极值性。但是对于许多物理或力学问题,要求场函数事先满足全部附加条件是很困难的。


参考: * Weak Solution - Wikipedia * 弹性力学变分提法随记 - 知乎

概述

我们可以从微分的概念出发理解变分。

微分:设有一个自变量 \(x\in(a,b)\) ,这个自变量的微小变化量记为 \(\mathrm{d}x\),称为微分。此时我们关心的是,当自变量有微小变化 (\(\mathrm{d}x\)) 时,函数 (\(y(x)\)) 变化了多少。

变分:设有一个函数(泛函) \(I=G(y,y^{\prime})\),其中 \(y=y(x)\) 是另一个函数(定义在某一个函数空间内是不确定的)。这个函数 \(y\) 的微小变化量记为 \(\delta y\),称为变分。此时我们关心的是,当函数 \(y\) 有微小变化时,函数(泛函) \(I\) 变化了多少。

简而言之,变分就是微分在函数空间的拓展,其精神内涵是一致的。

泛函:是将函数空间(无限维空间)映射到数域,就是把一个函数映射成一个数。打个比方,从A点到B点有无数条路径,每一条路径都是一个函数。这无数条路径,每一条函数(路径)的长度都是一个数。从这无数个路径当中选一个路径最短或者最长的,就是求泛函的极值问题。

变分

假设我们有两个定点 \((a,p)\)\((b,p)\),连接这两点的任意曲线的方程 \(y=y(x)\) 都将满足如下的边界条件:

\[ y(a)=p,\quad y(b)=q \]

现在考虑如下形式的定积分:

\[ I=\int_{a}^{b} f\left(y, y^{\prime}\right) \mathrm{d} x \]

其中 \(f(y, y^{\prime})\) 是关于 \(y(x)\) 和它的一阶导数 \(y^{\prime}(x)\) 的函数,在实际问题中,我们有时需要找到一个具体的 \(y(x)\) 使得 \(I\) 有极值。

注意在一般的极值问题中,我们考察的是自变量 \(x\) 的变化:\(x\) 取值多少时,函数会有极值。这个新问题的不同之处在于,我们考察的是函数 \(y(x)\) 的变化:\(y(x)\) 是什么形式时,\(I\) 会有极值(\(I\) 称作函数 \(y(x)\)泛函)。然而这两类问题依然有共通之处:当 \(I\) 取极值时,对 \(y(x)\) 作微小的变化,\(I\) 在一级近似下应该保持不变。

如果 \(y(x)\) 有微小改变 \(\delta y(x)\)(高大上叫法:\(\delta y(x)\) 称作函数 \(y(x)\) 的变分),那么 \(f\left(y, y^{\prime}\right)\) 的变化为:

\[ \delta f=\frac{\partial f}{\partial y} \delta y+\frac{\partial f}{\partial y^{\prime}} \delta y^{\prime} \]

\(I\) 相应的变化为:

\[ \delta I=\int_{a}^{b}\left[\frac{\partial f}{\partial y} \delta y+\frac{\partial f}{\partial y^{\prime}} \delta y^{\prime}\right] \mathrm{d} x \]

方括号里的第二项可以改写成 \(\frac{\partial f}{\partial y^{\prime}} \frac{\mathrm{d}(\delta y)}{\mathrm{d} x}\),然后可以进行分部积分:

\[ \int_{a}^{b} \frac{\partial f}{\partial y^{\prime}} \delta y^{\prime} \mathrm{d} x =\int_{a}^{b} \frac{\partial f}{\partial y^{\prime}} \mathrm{d}(\delta y) = \frac{\partial f}{\partial y^{\prime}} \delta y {\bigg|} _{a} ^{b} -\int _{a} ^{b} \delta y \frac{\mathrm{d}}{\mathrm{d} x}\left(\frac{\partial f}{\partial y^{\prime}}\right) \mathrm{d} x \]

由于 \(y(x)\) 的边界条件固定,\(\delta y(a)=\delta y(b)=0\),所以分部积分出来的第一项为零,仅第二项有贡献。上式化简为:

\[ \delta I=\int_{a}^{b}\left[\frac{\partial f}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{\partial f}{\partial y^{\prime}}\right)\right] \delta y(x) \mathrm{d}x \]

如果 \(I\) 有极值,对任意满足边界条件的 \(\delta y(x)\) 都必须有 \(\delta I=0\),这就要求:

\[ \frac{\partial f}{\partial y}-\frac{\mathrm{d}}{\mathrm{d} x}\left(\frac{\partial f}{\partial y^{\prime}}\right)=0 \]

这便是 Euler-Lagrange 方程,它是变分法的核心定理。有了它,原则上就可以找出所寻求的极值函数 \(y(x)\)

Euler-Lagrange 方程的应用:两点间最短路径问题、最速降曲线问题、悬链线等。


本文参考了以下博主分享的文章: * 浅谈变分原理 * 微分、差分和变分的概念有什么异同?

动词

引起学术界的广泛关注 generate a good deal of academic attention

名词

发生 occur 处理 process >> deal with > handle 大量关注 a lot of attention

形容词

各种各样的 various

副词

特别地 in particular > especially

连接词或短语

同时 at the same time >> meanwhile 此外 furthermore = moreover 在...的基础上 on the basis of 随后 in the fllowing

概述

主要区别: * 作用对象。FFT作用于时间序列信号;Bode图作用于系统传递函数 \(H(s)\) * 目的。FFT的目的是得到一段信号的频谱;Bode图的目的是得到传递函数的幅频特性和相频特性 * 原理。FFT基于傅里叶变换;Bode图基于拉普拉斯变换 * 适用性。FFT分析的时间序列可以来自于线性系统或非线性系统;Bode图只能分析线性非时变系统的传递函数

主要联系: * 两者都得到频域下的结果 * 在数学形式上,得到传递函数时使用的拉普拉斯变换退化后得到傅里叶变换

傅里叶变换

傅里叶变换的目的:得到被采集信号包含的频率成分和每个频率成分在原信号中的占比,用图形表示则为频谱图。

傅里叶变换简单通俗理解就是把看似杂乱无章的信号考虑成由一定振幅、相位、频率的基本正弦(余弦)信号组合而成,傅里叶变换的目的就是找出这些基本正弦(余弦)信号中振幅较大(能量较高)信号对应的频率,从而找出杂乱无章的信号中的主要振动频率特点。如减速机故障时,通过傅里叶变换做频谱分析,根据各级齿轮转速、齿数与杂音频谱中振幅大的对比,可以快速判断哪级齿轮损伤。

傅里叶正变换可以写为: \[ X(j\omega) = \int^{+\infty}_ {-\infty} x(t) e^{ -j \omega t} dt \]

但是,我们发现傅里叶变换中的基函数是 \(e^{j\omega t}\),能用它拟合的函数具有共同特征,那就是,各个频率下的三角函数分量不随时间衰减。这个假设具有明显的局限性,要么是理想没有阻尼的系统,要么是研究有阻尼的系统对特定频率的瞬时响应。

因此需要对傅里叶变换的基函数进行拓展,由此引出拉普拉斯变换:

\[ e^{\sigma t}e^{j\omega t}=e^{(\sigma+j\omega) t}=e^{st} \]

其中,\(s=\sigma+j\omega\) 为拉普拉斯变换中的变量,表示复实数(\(j\)为虚数单位)。

PS:有关傅里叶变换的详细内容可以参照另一篇博客:快速傅里叶变换(FFT)基本原理

拉普拉斯变换

拉普拉斯变换可以写为: \[ X(\sigma+j\omega)=\int^{+\infty}_ {-\infty} x(t) e^{(\sigma+j\omega) t} dt=\int^{+\infty}_ {-\infty} x(t) e^{s t} dt=X(s) \]

\[ X(s)=\int^{+\infty}_ {-\infty} x(t) e^{s t} dt \]

拓展傅里叶变换得到拉普拉斯变换。它是为简化计算而建立的实变量函数和复变量函数间的一种函数变换。对一个实变量函数作拉普拉斯变换,并在复数域中作各种运算,再将运算结果作拉普拉斯反变换来求得实数域中的相应结果,往往比直接在实数域中求出同样的结果在计算上容易得多。拉普拉斯变换的这种运算步骤对于求解线性微分方程尤为有效,它可把微分方程化为容易求解的代数方程来处理,从而使计算简化。在经典控制理论中,对控制系统的分析和综合,都是建立在拉普拉斯变换的基础上的。

引入拉普拉斯变换的一个主要优点,是可采用传递函数 \(H(s)\) 代替微分方程来描述系统的特性。这就为采用直观和简便的图解方法来确定控制系统的整个特性(见信号流程图、动态结构图)、分析控制系统的运动过程(见奈奎斯特稳定判据、根轨迹法),以及综合控制系统的校正装置(见控制系统校正方法)提供了可能性。

Bode图

刚才提到,通过拉普拉斯变换可以得到系统的传递函数;而传递函数可以直观的反应线性系统输入与输出之间的关系。当我们对传递函数可视化时往往就会用到Bode图

伯德图(bode)利用对数表示系统传递函数的幅频、相频特性,它分为两个图,一个是对数幅频图、一个是相频图,横坐标均为 \(\omega\) ,纵坐标一个为\(L(\omega)=20\lg A(\omega )\)表示对数幅值,一个为\(\varphi (\omega)\) 表示相位角。

:Bode图 假设系统的传递函数为: \[ H(s)=\frac{A}{s+a} \] \(s=\sigma+j\omega\),则传递函数可以写为: \[ H(\sigma+j\omega)=\frac{A}{(\sigma+j\omega)+a} \] 此时传递函数的图像是三维的,有两个自变量\(\sigma 和 \omega\),但是Bode图只关系其中的\(\omega\),将另一个变量视为0,\(\sigma=0\)\[ H(\omega)=\frac{A}{(j\omega)+a} \] 上式的对数幅值和相角分别为: \[ L(\omega)=20\lg \big\vert H(\omega) \big \vert = 20\lg \bigg( \frac{A}{\sqrt{\omega^2+a^2}} \bigg ) \] \[ \varphi (\omega) = - \arctan \frac{\omega}{a} \] 利用以上两式即可画出Bode图,横坐标均为\(\omega\)


本文参考了以下博主分享的文章: * 信号处理趣学D8——关于拉氏变换和频谱图的那些事儿 * 傅里叶变换、拉普拉斯变换、Z 变换的联系是什么?为什么要进行这些变换? * MATLAB频域分析,奈氏图、伯德图、对数幅相图绘制 * Bode plot - Wikipeida * 傅里叶变换和拉普拉斯变换的物理解释及区别 * Bode Plots(伯德图)

问题背景

在使用Matlab或Python解决某些实际问题时,我们不免要求解非线性方程组,或解决一些优化问题。这时你就有可能遇到 Levenberg–Marquardt 算法

Levenberg–Marquardt 算法

Levenberg–Marquardt 算法是一种最小二乘(模型拟合)算法。它可以解决以下问题:

\[ \min_{x}f(x)= \min_{x} \\| F(x) \\|^2_2 = \min_{x} \sum_i F_i^2(x) \]

即求一个 \(x\) 使 \(f(x)\) 最小。

1944年Kenneth Levenberg首次提出该算法,1963年Donald Marquardt也发现了该算法,后来这种算法被称作Levenberg–Marquardt。

关于该算法的详细原理和实现可以参考 Flannery 等人的著作 Numerical recipes in C

Levenberg–Marquardt 算法与非线性方程组求根

求解非线性方程组即求解一组 \(\bar{x}\) 令其满足 \(\bar{F}(x)=0\)。显然非线性方程组的根可以令上式成立。

因此,给定初始值,通过Levenberg–Marquardt 算法求得的 \(\bar{x}\) 就是非线性方程组 \(\bar{x}\) 的一组近似根(初始值附近)。

非线性方程组的根显然不止一组,而通过该方法只能得到距离初始值较近的根,要想得到全部的根,需要配合其他算法。

参考文献

  • Levenberg, K. (1944). A method for the solution of certain non-linear problems in least squares. Quarterly of applied mathematics, 2(2), 164-168.
  • Mardquardt, D. W. (1963). An algorithm for least square estimation of parameters. J. Soc. Ind. Appl. Math, 11, 431-441.
  • Flannery, B. P., Press, W. H., Teukolsky, S. A., & Vetterling, W. (1992). Numerical recipes in C. Press Syndicate of the University of Cambridge, New York, 24(78), 36.

问题背景

论文写作时,我们往往需要插入矢量图,以保证图形被阅览时不会因缩放而出现难看的锯齿(失真)。

需要矢量图的常用场景如下: * 纯矢量图 * 简单的折线、柱形、条形图等 * 抽象的模型图、概念图、流程图、程序框图等 * 混合类图形,混合了以上两种图形 * 矢量图与位图的混合 * 色彩复杂的模型图(位图)旁配上文字、简单的折线图、流程图等 * 在流程图中嵌入各种结果,如Science、Nature中常见的占用一页的超大混合类图形 * 其他

各种类型图片的区别

矢量图

简单来说矢量图就像用几何图形来描述一幅图,在矢量图放大时,我们所记录的几何图形的各种角度、形状等并没有改变,所以无论是放大还是缩小,都不会影响矢量图的清晰度。

文件类型: Adobe Illustrator的.AI.EPS.SVG、AutoCAD的.dwg.dxf、Corel DRAW的.cdr等。.pdf也常被用来保存矢量图。

制图: 使用Adobe Illustrator(AI)、CorelDRAW、FlashMX制图

位图

位图又被称作点阵图或栅格图像,它的特点就是,整幅图由许多的‘点’组成,这些‘点’我们称为“像素”。在位图模式下,计算机会将图片的每个像素点进行保存。当位图放大到一定程度时,我们会发现图片是由一个一个的小方块组成,这些小方块就是像素点。

文件类型: .jpg.png.bmp.pcx.gif.tif.psd等。

制图: 使用Photoshop、Painter和Windows系统自带的画图工具制图。

如何在论文中使用矢量图

这里主要讨论使用LaTex写作时,如何插入矢量图。

第一步:使用任何软件得到矢量图形。如果使用Matlab、Python得到折线、柱形图等,直接输出为.pdf.eps格式。如果使用Visio、PPT等软件绘制流程图、模型图,导出为.pdf。(貌似微软的办公软件不支持.eps

(如果是直接插入已得到的矢量图,请转至第三步)

第二步:使用Adobe Illustrator将多个矢量图组合在一起后用快捷键Shift+o裁剪画布,删去边框。将文件保存为.eps.pdf。如果只是删除.pdf中的白边,使用Adobe Acrobat 的编辑功能就能轻松做到。

第三步:在LaTex中插入组合好的矢量图。

提示:在第一步制作矢量图时,最好调整好字号,图宽以适合出现在论文中或PPT中,并且使用Matlab时一定要保存一份.fig文件,便于以后查找修改。

Office套件中插入矢量图

Office办公软件下插入矢量图,使用.emf即可。Office系列间的互通是较为方便的,如Visio、PPT、Word间相互插入文件。

使用Word和Latex一段时间以后,还是非常推荐大家使用LaTex的(强迫症福音)。毕竟插入的行内、行间公式怎是一个“服服帖帖”就能概括的。更别说切换模板时,只需要改动导言区就ok啦,简直不要太方便。LaTex缺点也有,表格制作很麻烦QAQ…,整体学习成本略高。

参考文献

  • https://www.zhihu.com/question/378251607/answer/1160351444

引言

快速傅里叶变换(FFT)是离散傅里叶变换(DFT)的加速算法,而DFT则是将连续的傅里叶变换离散化(在时域和频域离散),连续傅里叶变换可由傅里叶展开式推导得出。

正向逻辑为: 傅里叶展开式 \(\rightarrow\) 推导得到 \(\rightarrow\) 连续傅里叶变换对 \(\rightarrow\) 离散后得到 \(\rightarrow\) 离散傅里叶变换 \(\rightarrow\) 开发加速算法 \(\rightarrow\) 快速傅里叶变换 ## 傅里叶变换 ### 三角函数形式 对于一个周期为 \(T\) 的信号 \(x(t)\),其基频为 \(f_1=1/T\)。对 \(x(t)\) 进行傅里叶展开:

\[ x(t) = \frac{a_0}{2} +\sum^{\infty}_ {n=1} \left[ a_n \cos\left(2\pi n f_1 t \right) + b_n \sin\left(2\pi n f_1 t \right) \right] \]

上式中的 \(a_n\)\(b_n\) 是傅里叶系数:

\[ a_n = \frac{2}{T} \int^T_ 0 x(t) \cos \left( 2\pi n f_1 t \right) dt \]

\[ b_n = \frac{2}{T} \int^T_ 0 x(t) \sin \left( 2\pi n f_1 t \right) dt \]

指数形式

根据欧拉公式,可以将傅里叶变换式写为指数形式,指数形式看起来更为简洁,利于推导。 根据欧拉公式,有以下变换:

\[ \cos \left( 2\pi n f_1 t \right) = \frac{1}{2} \left[ e^{j2\pi n f_1 t} + e^{-j2\pi n f_1 t} \right] \]

\[ \sin \left( 2\pi n f_1 t \right) = \frac{1}{2j} \left[ e^{j2\pi n f_1 t} - e^{-j2\pi n f_1 t} \right] \]

欧拉公式:\(e^{ix} = \cos x + i \sin x\)

带入傅里叶展开式,可得:

\[ x(t) = \frac{a_0}{2} +\sum^{\infty}_ {n=1} \left[ \frac{a_n}{2} \left[ e^{ j 2\pi n f_1 t } + e^{ -j 2\pi n f_1 t } \right] + \frac{b_n}{2j} \left[ e^{ j 2\pi n f_1 t } - e^{ -j 2\pi n f_1 t } \right] \right] \]

合并同类项可得:

\[ x(t) = \frac{a_0}{2} +\sum^{\infty}_ {n=1} \left[ \frac{a_n - j b_n}{2} e^{ j 2\pi n f_1 t } + \frac{a_n + b_n}{2} e^{ -j 2\pi n f_1 t } \right] \]

观察指数项前的系数,令

\[ X(2\pi n f_1) = \frac{a_n - j b_n}{2} \]

则有:

\[ X(-2\pi n f_1) = \frac{a_n + j b_n}{2} \]

发现指数项前的系数具有奇偶性,利用奇偶性并改变求和符号的下限,将展开式改写为:

\[ x(t) = \sum^{\infty}_ {n=-\infty} X(2\pi n f_1) e^{ j 2\pi n f_1 t } \]

注意,此时求和符号的下限由一变为负无穷,\(1 \rightarrow -\infty\)。式中的 $X(2n f_1) $ 展开后为:

\[ \begin{gathered} X(2\pi n f_1) =\frac{1}{2} \left[ \frac{2}{T} \int^T_ 0 x(t) \cos \left( 2\pi n f_1 t \right) dt - \frac{2j}{T} \int^T_ 0 x(t) \sin \left( 2\pi n f_1 t \right) dt \right] \\\\ = \frac{1}{T} \int^T_ 0 x(t) \left[ \cos \left( 2\pi n f_1 t \right) - j \sin \left( 2\pi n f_1 t \right) \right] dt \end{gathered} \]

再利用欧拉公式,可得:

\[ X(2\pi n f_1) = \frac{1}{T} \int^T_ 0 x(t) e^{ -j 2\pi n f_1 t} dt \]

总结上述推导过程,得到指数形式的傅里叶变换对:

\[ x(t) = \sum^{\infty}_ {n=-\infty} X(2\pi n f_1) e^{ j 2\pi n f_1 t } \]

\[ X(2\pi n f_1) = \frac{1}{T} \int^T_ 0 x(t) e^{ -j 2\pi n f_1 t} dt \]

以上两个公式中,第二式才是傅里叶变换,第一式是傅里叶变换反演公式。

当信号没有周期,即周期无穷大时 \(T\rightarrow \infty\),有\(f_1\rightarrow 0\),此时第一式的求和符号可以改写为积分号,\(f_1\) 可以被当作是一个连续变化的量,傅里叶变换对就被改写为了连续形式。

傅里叶变换的实际意义

傅里叶变换对第一式:表示被采集信号 \(x(t)\) 可以由具有不同频率(\(2\pi n f_1\))和幅值(\(X(2\pi n f_1)\))的周期信号叠加得到。

傅里叶变换对第二式:为傅里叶变换,表示不同频率分量前的系数,可以理解为,该频率分量在被采集信号中的占比。

傅里叶变换的目的:得到被采集信号包含的频率成分和每个频率成分在原信号中的占比,用图形表示则为频谱图。

离散傅里叶变换

离散傅里叶变换对应的是在时域、频域都有限长,且都是离散的一类变换。

现有一个时域模拟信号 \(x(t)\),对其进行采样,采样长度为 \(T\), 采样点有 \(N\) 个,采样时间间隔为 \(\Delta t\)

则采样长度与采样点的关系为:

\[ T=N\Delta t \]

采样频率为:

\[ f_s = \frac{1}{\Delta t} \]

频谱的基频是\(f_1\),即:

\[ f_1=\frac{1}{T} \]

根据傅里叶变换公式,由傅里叶变换求得的频谱谱线都是基频的整数倍,即频谱谱线的间隔:

\[ \Delta f=f_1=\frac{1}{T} = \frac{1}{N\Delta t} \]

记采样得到的离散的时域数字信号为\(x(k\Delta t)\),把傅里叶变换式中的连续变量替换为离散变量

\[ t=k\Delta t,\quad 2\pi nf_1=2\pi n \Delta f \]

再将\(T=N\Delta t\)带入傅里叶变换对,对时间积分变为有限时间段内求和,即,\(\int \rightarrow \sum\)\(dt\rightarrow\Delta t\),得到:

\[ x(k\Delta t) = x(k)= \sum^{N-1}_ {n=0} X(2\pi n \Delta f) e^{ j 2\pi n k / N } \]

\[ X(2\pi n \Delta f) =X(n)= \frac{1}{N} \sum^{N-1}_ {k=0} x(k\Delta t) e^{ -j 2\pi n k /N} \]

习惯上将标定因子 \(N\) 移至反变换式中去,并且用 \(n\) 表示第 \(n\) 个采样点,用 \(k\) 表示第 \(k\) 条谱线(频率分量),即将上式中 \(n\)\(k\) 的位置互换。总结上述结论有:

\[ x(n) = \frac{1}{N} \sum^{N-1}_ {n=0} X(k) e^{ j 2\pi n k / N } \]

\[ X(k) = \sum^{N-1}_ {k=0} x(n) e^{ -j 2\pi n k /N} \]

上式就是离散傅里叶变换式。式中:

\[ k= 0,1,2\dots,N-1; \quad n= 0,1,2\dots,N-1 \]

根据离散傅里叶变换式可以求出第 \(n\) 条谱线对应的傅里叶系数的值,即该频率分量在整个信号中的占比。

我们注意到离散傅里叶变换式中求和符号的上下限改变了,不再是正负无穷。再改为正负无穷可以吗?改了之后两个公式含义还一样吗?答案是一样的,要用采样定理回答该问题。虽然标注是无穷,但是采样频率一定要大于2倍的被采样信号最高频率:

\[ f_s=\frac{1}{\Delta t}>2f_m \]

当采样频率已经确定为 \(1/ \Delta t\) 时,满足采样定理的 \(x(t)\) 的最高频率分量是:

\[ f_m<\frac{1}{2}f_s=\frac{1}{2\Delta t} = \frac{1}{2} N \Delta f \]

即满足采样定理的最后一根谱线所在的频率为\(1/(2N\Delta f)\),换而言之 \(k\) 最大为 \(N/2\)

快速傅里叶变换

FFT提出的背景

使用离散傅里叶变换时,计算一点的 \(X(n)\) 需要 \(N\) 次复数乘法和 \(N-1\) 次复数加法运算;计算全部的频谱需要 \(N^2\) 次复数乘法和 \(N(N-1)\) 次复数加法。而实现1次复数乘法需要4次实数乘法和2次实数加法,实现1次复数加法,需要2次实数加法。\(N=1024\)时,需要 \(1048576\) 次复数乘法运算。

有没有什么方法可以缩短计算时间呢?

1965年J. W. Cooley 和 J. W. Tukey 提出了一种快速求解DFT的算法,将乘法运算量由 \(N^2\) 降低至 \(\frac{N}{2} \log_ 2 N\) 次。以 \(N=1024\) 为例,他们提出的算法,只需要 \(5120\) 次复数乘法运算。

值得注意的是FFT不是特指某一种算法,而是指一类算法,自1965年后也有许多优秀学者提出了新的FFT算法。

FFT基本原理

对于 \(N\) 点序列 \(x(n)\),其 DFT变换对为:

\[ x(n) = \frac{1}{N} \sum^{N-1}_ {n=0} X(k) W^{-nk}_ N \]

\[ X(k) = \sum^{N-1}_ {k=0} x(n) W^{nk}_ N \]

其中,\(W_ N = e^{-j2\pi /N}\)

DFT还可以写成矩阵形式,定义 \(\mathbf{W}_N\) 为:

\[ \mathbf{W}_N = \left[ W^{nk} \right]= \begin{bmatrix} W^0 & W^0 & W^0 & \cdots & W^0 \\\\ W^0 & W^1 & W^2 & \cdots & W^{N-1} \\\\ W^0 & W^2 & W^4 & \cdots & W^{2(N-1)} \\\\ \vdots & \vdots & \vdots & \vdots & \vdots \\\\ W^0 & W^{(N-1)} & W^{2(N-1)} & \cdots & W^{(N-1)(N-1)} \end{bmatrix} \]

\[ \bar{X}_N = \begin{bmatrix} X(0) & X(1) & \cdots & X(N-1) \end{bmatrix}^\mathrm{T} \]

\[ \bar{x}_N = \begin{bmatrix} x(0) & x(1) & \cdots & x(N-1) \end{bmatrix}^\mathrm{T} \]

则DFT的矩阵表达式为:

\[ \bar{X}_N = \mathbf{W}_N \bar{x}_N \]

观察矩阵形式的DFT表达式,可以发现DFT中包含大量重复运算,矩阵 \(\mathbf{W}_N\) 虽然有 \(N^2\) 个元素,但是其中只有 \(N\) 个独立值,并且一部分元素的值极为简单。\(W\) 因子的取值有如下特点:

  • \(W^0=1\), \(W^{N/2}=-1\)
  • \(W^{N+r}=W^r\), \(W^{N/2+r}=-W^r\)

对于4点的DFT,需要 \(4^2=16\) 次复数乘法,但是利用 \(W\) 因子的周期性和对称性可以大大简化计算,相关的DFT矩阵格式为:

\[ \begin{bmatrix} X(0) \\\\ X(1) \\\\ X(2) \\\\ X(3) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 \\\\ 1 & W^1 & -1 & -W^1 \\\\ 1 & -1 & 1 & -1 \\\\ 1 & -W^1 & -1 & W^1 \end{bmatrix}\begin{bmatrix} x(0) \\\\ x(1) \\\\ x(2) \\\\ x(3) \end{bmatrix} \]

将矩阵的第二列和第三列交换,可得:

\[ \begin{bmatrix} X(0) \\\\ X(1) \\\\ X(2) \\\\ X(3) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 \\\\ 1 & -1 & W^1 & -W^1 \\\\ 1 & 1 & -1 & -1 \\\\ 1 & -1 & -W^1 & W^1 \end{bmatrix}\begin{bmatrix} x(0) \\\\ x(2) \\\\ x(1) \\\\ x(3) \end{bmatrix} \]

由此可得:

\[ \begin{gathered} X(0) = [x(0) + x(2)] + [x(1) + x(3)] \\\\ X(1) = [x(0) - x(2)] + [x(1) - x(3)] W^1 \\\\ X(2) = [x(0) + x(2)] - [x(1) + x(3)] W^1 \\\\ X(3) = [x(0) - x(2)] - [x(1) - x(3)] \end{gathered} \]

利用上式计算DFT只需要 \(1\) 次复数乘法运算。处理长序列时,只需要将长序列分成类似于上式或比上式更简短的序列后,进行简单的运算,再按一定方式组合起来即可。

因此,FFT的基本原理是:先将原始的序列分解为一系列的短序列,充分利用 \(W\) 因子的周期性和对称性,进而求出这些短序列相应的DFT并进行适当组合,达到删除重复计算,减少乘法运算和简化结构的目的。

对于大部分FFT算法,有一些通用的概念和规律: * :将原信号每次折半,分为更小的单元,每折半一次,多出一级 * 蝶形单元:在FFT计算过程中,计算并不是顺序的,画出信号流图会发现,信号流图中包含着大量8字型(类似于蝴蝶)的计算单元 * :每一级的蝶形单元,按照其特性可以分为若干组 * \(W\)因子的分布:每一级的 \(W\) 因子具有一定的分布规律 * 码位倒置:使用FFT时,输出序列 \(\bar{X}_N\) 依照正序排列,但是输入序列的次序不再是自然序列,其排布次序和二进制码翻转、二进制与十进制转换有关

本文不再详细叙述具体的FFT算法。

0%