CN117350112A - Fokker-Planck方程的无限元求解方法 - Google Patents
Fokker-Planck方程的无限元求解方法 Download PDFInfo
- Publication number
- CN117350112A CN117350112A CN202311289338.9A CN202311289338A CN117350112A CN 117350112 A CN117350112 A CN 117350112A CN 202311289338 A CN202311289338 A CN 202311289338A CN 117350112 A CN117350112 A CN 117350112A
- Authority
- CN
- China
- Prior art keywords
- probability density
- function
- equation
- density function
- fokker
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 85
- 239000011159 matrix material Substances 0.000 claims abstract description 62
- 238000004364 calculation method Methods 0.000 claims abstract description 21
- 230000014509 gene expression Effects 0.000 claims abstract description 21
- 238000009792 diffusion process Methods 0.000 claims abstract description 16
- 238000006073 displacement reaction Methods 0.000 claims abstract description 15
- 238000010276 construction Methods 0.000 claims abstract description 9
- 239000000654 additive Substances 0.000 claims abstract description 8
- 230000000996 additive effect Effects 0.000 claims abstract description 8
- 230000005284 excitation Effects 0.000 claims description 12
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 230000004044 response Effects 0.000 claims description 4
- 102100029469 WD repeat and HMG-box DNA-binding protein 1 Human genes 0.000 claims description 2
- 101710097421 WD repeat and HMG-box DNA-binding protein 1 Proteins 0.000 claims description 2
- 230000002146 bilateral effect Effects 0.000 claims description 2
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 239000006185 dispersion Substances 0.000 claims description 2
- 230000004907 flux Effects 0.000 claims description 2
- 238000010606 normalization Methods 0.000 claims description 2
- 230000003595 spectral effect Effects 0.000 claims description 2
- 230000008878 coupling Effects 0.000 claims 1
- 238000010168 coupling process Methods 0.000 claims 1
- 238000005859 coupling reaction Methods 0.000 claims 1
- 238000012163 sequencing technique Methods 0.000 claims 1
- 238000004458 analytical method Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000013016 damping Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 239000000725 suspension Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2101/00—Indexing scheme relating to the type of digital function generated
- G06F2101/14—Probability distribution functions
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Operations Research (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Computing Systems (AREA)
- Complex Calculations (AREA)
Abstract
Fokker‑Planck方程的无限元求解方法，包括：由加性高斯白噪声激励的机械动力学系统运动控制方程，推导出含概率密度函数的Fokker‑Planck方程；无限元方法中，几何插值函数使用双线性插值函数，场插值函数使用衰减型无限元插值函数，推导单元系数矩阵m和k的表达式，根据有限元矩阵构造原理得到全局系数矩阵M和K；使用算子分裂法对单元系数矩阵k进行拆分，将其分裂为对流项和扩散项，得到单元系数矩阵k对流项和扩散项的具体表达式；给定参数下，绘出概率密度函数随时间演化的等高线图和稳态时刻概率密度函数沿位移、速度方向的截面图。本发明弥补了有限元方法中边界处因人为截断而丢失的概率密度函数值，在一定计算精度下，结合算子分裂法使得计算时间大大降低。
Description
技术领域
本发明属于随机动力学系统求解方法和技术领域，具体涉及一种无限元法结合算子分裂法求解Fokker-Planck(简称FP，中文翻译为福克-普朗克)方程的数值方法。
背景技术
自然、工程及社会领域中存在许多随机振动现象，诸如地震时大地的随机振动、大气湍流对飞机的随机扰动、海浪对轮船的随机冲击、路面不平度的随机分布激励对车辆的随机振动等。这些随机振动都不能用时间与空间坐标的确定性函数来描述，在动力学响应分析中通常用概率或统计的方法进行描述和建模。Fokker-Planck方程控制一类特定机械动力学系统(线性或非线性系统)在随机激励作用下的动态响应的概率密度函数，包含系统在状态空间中任何时刻的完整概率密度分布。可见，Fokker-Planck方程在随机振动中发挥着至关重要的作用。
然而存在精确解的Fokker-Planck方程很少，只有在严格限制条件下，对极少数情况才能求得精确解，因此数学家们研究了很多数值方法求解Fokker-Planck方程，其中有限元法最为准确。1981年Bergman首次提出借助Petrov-Galerkin有限元法求解后向Fokker-Planck方程，解决了在加性高斯白噪声激励下二阶线性系统的首次穿越问题。2013年Pichler等研究了有限元法和有限差分法的区别，通过数值仿真算例验证了有限元法的求解精度更高。但是随着机械动力学系统自由度的增加，计算维度不断增长，待求解的标量方程数也会呈现指数增长，使用传统的有限单元法求解Fokker-Planck方程需要花费大量的计算时间。
为了快速、准确的求解Fokker-Planck方程，本发明提出了无限元法结合算子分裂方法，在传统的有限域中耦合无限域，弥补了有限元法求解边界处因人为截断而丢失的概率密度函数值，结合算子分裂法对每个子部分进行求解，在一定计算精度下使得计算时间为有限单元法的1～5％。
发明内容
本发明在于求解由加性高斯白噪声激励的机械动力学系统的Fokker-Planck方程。目的在于以无限元法结合算子分裂法，在一定计算精度下，快速求解机械动力学系统的Fokker-Planck方程，获得机械动力学系统任意时刻的完整概率密度函数值。
实现本发明目的的技术解决方案为：Fokker-Planck方程的无限元求解方法，包括以下步骤：
步骤1：由加性高斯白噪声激励的机械动力学系统运动控制方程，推导出含概率密度函数的Fokker-Planck方程；
步骤2：无限元方法中，几何插值函数使用双线性插值函数，场插值函数使用衰减型无限元插值函数。推导单元系数矩阵m和k的表达式，根据有限元矩阵构造原理得到全局系数矩阵M和K；
步骤3：使用算子分裂法对单元系数矩阵k进行拆分，将其分裂为对流项和扩散项，得到单元系数矩阵k对流项和扩散项的具体表达式；
步骤4：给定参数下，绘出概率密度函数随时间演化的等高线图和稳态时刻概率密度函数沿位移、速度方向的截面图。最终无限元方法求解Fokker-Planck方程得到的解是一个维度为10201的矩阵，矩阵中全部的概率密度函数值分布在由状态变量位移和速度组成的状态空间中，且稳态时刻的概率密度函数值收敛于高斯型分布。
本发明以无限元法为基础，弥补了有限元方法中边界处因人为截断而丢失的概率密度函数值，在一定计算精度下，结合算子分裂法使得计算时间大大降低，为Fokker-Planck方程的求解提供新的高效求解方法。
本发明与现有技术相比，其显著优点为：
(1)考虑了有限元方法中边界处因人为截断而丢失的概率密度函数值，使得计算结果更精确。
(2)使用算子分裂法对单元系数矩阵k进行拆分时，将其按照对流过程和扩散过程拆分为对流项和扩散项，在一定计算精度下，花费更少的计算时间得到系统任意时刻概率密度函数值。
附图说明
图1为本方法实施步骤流程图；
图2为无限元网格图；
图3为稳态时概率密度尖端点沿位移方向横截面图；
图4为稳态时概率密度尖端点沿速度方向横截面图；
图5为原点概率密度函数随归一化时间的演化图；
图6为概率密度函数随时间的演化等高线图；
图7为概率密度函数稳态时刻三维图。
具体实施方式
本发明在无限元方法基础上结合算子分裂方法，在一定计算精度下大大减少了求解机械动力学系统Fokker-Planck方程所需的时间。下面将结合附图对本发明专利的技术方案进行完整地描述，包括以下步骤：
步骤1：由加性高斯白噪声激励的机械动力学系统运动控制方程，推导出含概率密度函数的Fokker-Planck方程；
由高斯白噪声v(t)激励的机械动力学系统g(x)的控制方程可写成如下形式
其中x是二维域实数域R2中动力学系统的响应，通过定义向量h可以控制机械动力学系统的振动类型。在方程(1)右侧存在的高斯白噪声激励v(t)由其前两阶导数矩控制，一阶矩表示均值，二阶矩表示协方差函数，如方程：
E[v(t)]＝0 (2)
E[v(t)v(t+t′)]＝2Dδ(t′) (3)
E[·]是均值函数，δ(·)是狄拉克函数，D/π表示输入的恒定双边谱密度函数的幅值。
在高斯白噪声激励的所有统计量中，其高于二阶的导数矩都为0，因此系统的特征完全由激励的前两阶矩控制，可以得到由第一阶矩和第二阶矩定义的Fokker-Planck方程：
其中概率密度函数p(x,t)在二维实数域中满足归一化条件；
在求解Fokker-Planck方程时需要定义初始条件和边界条件，方程的解通常是计算域中的局部分布，大部分边缘处的概率密度函数的值为0，因此在位移方向选择零通量边界条件。初始条件使用双正态概率密度函数：
其中φ(·)是正态概率密度函数，μx01和σ01表示x1方向的期望和标准差，μx02和σ02表示x2方向的期望和标准差。
步骤2：无限元方法中，几何插值函数使用双线性插值函数，场插值函数使用衰减型无限元插值函数。推导单元系数矩阵m和k的表达式，根据有限元矩阵构造原理得到全局系数矩阵M和K；
推导方程的弱解形式，将方程(4)右侧所有项移到方程的左侧，每项都乘以与概率密度p同型的权函数w，然后方程两侧同时在计算域Ω上积分，为了表达更直观，将g1(x)、g2(x)简写为g1、g2；
积分计算域Ω为±∞，公式(7)可以利用分部积分进行化简，再结合几何边界条件，边界处的概率密度的值为0，可以得到机械动力学系统控制方程的弱解形式：
在有限单元中，使用双线性插值函数对概率密度离散。有限元计算域Ω被分为n个小单元，每个单元具有四个节点，如图2中红色区域所示。每个节点具有一个自由度，单元内部均由双线性插值计算，并通过坐标转换公式转为等参单元进行计算。因为权重函数w(x,t)和概率密度函数p(x,t)形式相同，所以可以用插值函数表示为：
公式(9)和公式(10)中，Nr(x)、Ns(x)为权函数w和概率密度函数p在某个单元对应的双线性插值函数，和/>为节点处的概率密度函数值，r,s＝1,2,3,4。给出双线性插值函数的具体形式，其中-1≤ξ≤1,-1≤η≤1
将公式(9)和公式(10)带入公式(8)，化简可得有限元离散系统的单元运动方程：
其中单元系数矩阵m和k的具体形式为：
m＝∫∫NNTdx1dx2 (16)
N表示全局插值函数矩阵。
得到所有的单元系数矩阵后，根据节点编号和单元编号对应的方式构造全局系数矩阵M和K，从而全局运动控制方程为：
图2中黑色区域为x2正向和负向扩展到无穷远的区域，这两部分使用无限元方法进行分析，几何插值函数仍使用双线性插值函数，场插值函数使用衰减型插值函数进行插值离散。x2正向场插值函数为
公式(19)～(22)中，-1≤ξ≤1,0≤η≤∞，ηi和L为决定衰减程度的任意常数，可以改变ηi/L的值来设定不同的衰减程度。
x2负方向场插值函数为：
公式(23)～(26)中，-1≤ξ≤1,-∞≤η≤0。
无限元方法中场插值函数不同与无限元方法中的插值函数，因此权重函数和概率密度函数的离散放表达式也不同，x2正向表示为：
Nr′(x)、Ns′(x)为x2正向权函数w和概率密度函数p在某个单元对应的衰减型插值函数。
x2负向表示为：
Nr″(x)、Ns″(x)为x2负向权函数w和概率密度函数p在某个单元对应的衰减型插值函数。
无限域中的单元系数矩阵的表达式与公式(16)、(17)相似，只需把有限元法中的全局插值函数矩阵N换为无限元方法中的全局插值矩阵N′和N″。将有限域和无限域耦合，构造全局系数矩阵时同样按照节点编号和单元编号一致的原则，按照从左到右，从上到下的原则进行排序。为了构造时逻辑清晰，所以先对正向无限域进行组装，然后是有限元域，最后是负向无限域，由上到下获得全局系数矩阵M和K。
步骤3：使用算子分裂法对单元系数矩阵k进行拆分，将其分裂为对流项和扩散项，得到单元系数矩阵k对流项和扩散项的具体表达式；
借助无限元方法，得到每个部分的单元系数矩阵，然后结合算子分裂方法，将单元系数k分裂为k1和k2两项。公式(15)可以写成
其中k1与对流项相关，k2与扩散项相关。k1和k2具体的表达式如下
根据同样的矩阵构原理，得到全局系数矩阵K1和K2，全局运动控制方程可以写为：
步骤4：给定参数下，绘出概率密度函数随时间演化的等高线图和稳态时刻概率密度函数沿位移、速度方向的截面图。最终无限元方法求解Fokker-Planck方程得到的解是一个维度为10201的矩阵，矩阵中全部的概率密度函数值分布在由状态变量位移和速度组成的状态空间中，且稳态时刻的概率密度函数值收敛于高斯型分布。
使用算子分裂方法求解时，在时间区域[0,T]内，将其划分为n个Δt时间步长。在第i个时间步长内，使用算子分裂连续求解如下两个方程
在求解方程(35)时，以q0＝pi作为初值计算一个时间步长Δt，求得一个解qi+1；在求方程(36)时，以p0＝qi+1作为初值计算一个时间步长Δt，求得ti+1时刻的概率密度值pi+1。按照以上求解步骤循环求解，最终可以求得每个时刻的概率密度函数值。
实施例1
工程领域中存在诸多随机振动现象，例如汽车在路面上行驶时，路面不平度的随机分布激励引起车辆的随机振动、飞机降落时，起落架受到地面的随机激励从而引起飞机的随机振动，可以用机械动力学系统中的线性系统完成对四分之一车模型、飞机起落架进行分析和建模，同时这些随机振动可以使用线性系统的Fokker-Planck进行描述，然后借助无限元法结合算子分裂法求解系统的Fokker-Planck方程，从而实现机械动力学系统随机振动的分析。
下面结合具体实施案例对本发明进行进一步的介绍，应用本发明提出的一种Fokker-Planck方程的无限元求解方法的单自由度线性系统随机振动的分析方法，具体方法如下：
步骤1：由加性高斯白噪声激励的机械动力学系统运动控制方程，推导出含概率密度函数的Fokker-Planck方程；
对高斯白噪声激励下单自由度弹簧阻尼线性系统进行分析，给出系统的运动方程
其中ζ表示系统阻尼比，ω0表示系统的固有频率，x(t)表示状态空间中的随机过程。
令x(t)等于x1，等于x2，化简后带入方程(7)中得到机械动力学系统的Fokker-Planck方程
单自由度线性系统参数和Fokker-Planck方程的初始条件参数如表1所示
表1线性系统参数和初始条件参数
步骤2：无限元方法中，几何插值函数使用双线性插值函数，场插值函数使用衰减型无限元插值函数。推导单元系数矩阵m和k的表达式，根据有限元矩阵构造原理得到全局系数矩阵M和K；
步骤3：使用算子分裂法对单元系数矩阵k进行拆分，将其分裂为对流项和扩散项，得到单元系数矩阵k对流项和扩散项的具体表达式；
步骤4：给定参数下，绘出概率密度函数随时间演化的等高线图和稳态时刻概率密度函数沿位移、速度方向的截面图。最终无限元方法求解Fokker-Planck方程得到的解是一个维度为10201的矩阵，矩阵中全部的概率密度函数值分布在由状态变量位移和速度组成的状态空间中，且稳态时刻的概率密度函数值收敛于高斯型分布。
给定时间序列t为[0,20]，步长Δt＝0.01，x1的范围是(-10,10)，x2的范围是(-10,10)，衰减常数L＝5,ηi＝25，图3和图4为稳态时刻概率密度函数沿位移、速度方向的截面图，图6所示为线性系统概率密度函数随时间演化的等高线图。
步骤5：步骤4中的图6是单自由度系统的概率密度函数等高线图，几个归一化时间步长的叠加来显示单自由度系统的随机振动演化过程。从图中每个时刻的等高线可以看出每个时间步长的等高线并不均匀，说明系统的随机振动具有不规律性，越密集的等高线表明概率密度函数越集中，表明此时刻系统的振动越强烈，当归一化时间在1附近时，图像完成一次翻转，此时单自由度系统完成一个方向的振动。此外等高线图中每半个周期的概率密度函数具有对称性，说明单自由度系统在受随机激励后的振动沿着两个方向往复进行，最终的概率密度函数趋于稳态，此时系统的振动停止，达到稳定状态。
实施例2
本实施例提供应用本发明的无限元法结合算子分裂法求解Fokker-Planck方程的数值计算方法的机械动力学系统非线性振动评估方法，具体包括：
步骤1：由加性高斯白噪声激励的机械动力学系统运动控制方程，推导出含概率密度函数的Fokker-Planck方程；
步骤2：无限元方法中，几何插值函数使用双线性插值函数，场插值函数使用衰减型无限元插值函数。推导单元系数矩阵m和k的表达式，根据有限元矩阵构造原理得到全局系数矩阵M和K；
步骤3：使用算子分裂法对单元系数矩阵k进行拆分，将其分裂为对流项和扩散项，得到单元系数矩阵k对流项和扩散项的具体表达式；
步骤4：给定参数下，绘出概率密度函数随时间演化的等高线图和稳态时刻概率密度函数沿位移、速度方向的截面图。最终无限元方法求解Fokker-Planck方程得到的解是一个维度为10201的矩阵，矩阵中全部的概率密度函数值分布在由状态变量位移和速度组成的状态空间中，且稳态时刻的概率密度函数值收敛于高斯型分布。
步骤5：根据机械动力学系统的概率密度函数，评估机械动力学系统受随机激励后的随机振动。从图6机械动力学系统的概率密度函数等高线图可以看出，不同时刻的等高线密集程度不同，在系统振动的初始时刻较集中，可以评估机械动力学系统受随机激励后在初始时刻的随机振动较强烈；图5是概率密度函数随归一化时间的演化图，在整个时间区间内，概率密度函数一致存在且时刻发生变化，说明机械动力学系统受激励后一直存在振动，从图7概率密度函数的稳态图，从我们可以得出结论，此时机械动力学系统的非线性振动达到稳态。随着现代机械结构的日益复杂，机械动力学系统的非线性振动也难以预测，应用本发明的无限元法结合算子分裂法求解Fokker-Planck方程的数值计算方法的机械动力学系统非线性振动评估方法，可以对车辆悬架系统、飞机起落架等机械动力学系统的非线性振动进行评估。
本实施例的步骤1-4包含本发明技术方案的步骤1-4的所有技术手段。
本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举。本发明的保护范围不应当被视为仅限于实施例所陈述的具体结构和参数，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。
Claims (5)
1.Fokker-Planck方程的无限元求解方法，其特征在于，包括如下步骤：
步骤1：由加性高斯白噪声激励的机械动力学系统运动控制方程，推导出含概率密度函数的Fokker-Planck方程；
步骤2：无限元方法中，几何插值函数使用双线性插值函数，场插值函数使用衰减型无限元插值函数。推导单元系数矩阵m和k的表达式，根据有限元矩阵构造原理得到全局系数矩阵M和K；
步骤3：使用算子分裂法对单元系数矩阵k进行拆分，将其分裂为对流项和扩散项，得到单元系数矩阵k对流项和扩散项的具体表达式；
步骤4：给定参数下，绘出概率密度函数随时间演化的等高线图和稳态时刻概率密度函数沿位移、速度方向的截面图。最终无限元方法求解Fokker-Planck方程得到的解是一个维度为10201的矩阵，矩阵中全部的概率密度函数值分布在由状态变量位移和速度组成的状态空间中，且稳态时刻的概率密度函数值收敛于高斯型分布。
2.如权利要求1所述的Fokker-Planck方程的无限元求解方法，其特征在于，步骤1具体包括：
由高斯白噪声v(t)激励的机械动力学系统g(x)的控制方程可写成如下形式
其中x是二维域实数域R2中动力学系统的响应，通过定义向量h可以控制机械动力学系统的振动类型；在方程(1)右侧存在的高斯白噪声激励v(t)由其前两阶导数矩控制，一阶矩表示均值，二阶矩表示协方差函数，如方程
E[v(t)]＝0 (2)
E[v(t)v(t+t′)]＝2Dδ(t′) (3)
E[·]是均值函数，δ(·)是狄拉克函数，D/π表示输入的恒定双边谱密度函数的幅值；
在高斯白噪声激励的所有统计量中，其高于二阶的导数矩都为0，因此系统的特征完全由激励的前两阶矩控制，可以得到由第一阶矩和第二阶矩定义的Fokker-Planck方程
其中概率密度函数p(x,t)在二维实数域中满足归一化条件
在求解Fokker-Planck方程时需要定义初始条件和边界条件，方程的解通常是计算域中的局部分布，大部分边缘处的概率密度函数的值为0，因此在位移方向选择零通量边界条件；初始条件使用双正态概率密度函数：
其中φ(·)是正态概率密度函数，μx01和σ01表示x1方向的期望和标准差，μx02和σ02表示x2方向的期望和标准差。
3.如权利要求1所述的Fokker-Planck方程的无限元求解方法，其特征在于，步骤2具体包括：
推导方程的弱解形式，将公式(4)右侧所有项移到方程的左侧，每项都乘以与概率密度p同型的权函数w，然后方程两侧同时在计算域Ω上积分，为了表达更直观，将g1(x)、g2(x)简写为g1、g2
积分计算域Ω为±∞，公式(7)可以利用分部积分进行化简，再结合几何边界条件，边界处的概率密度的值为0，可以得到机械动力学系统控制方程的弱解形式
在有限单元中，使用双线性插值函数对概率密度离散；有限元计算域Ω被分为n个小单元，每个单元具有四个节点，每个节点具有一个自由度，单元内部均由双线性插值计算，并通过坐标转换公式转为等参单元进行计算；因为权重函数w(x,t)和概率密度函数p(x,t)形式相同，所以可以用插值函数表示为
公式(9)和公式(10)中，Nr(x)、Ns(x)为权函数w和概率密度函数p在某个单元对应的双线性插值函数，和/>为节点处的概率密度函数值，r,s＝1,2,3,4；给出双线性插值函数的具体形式，其中-1≤ξ≤1,-1≤η≤1
将公式(9)和公式(10)带入公式(8)，化简可得有限元离散系统的单元运动方程：
其中单元系数矩阵m和k的具体形式为
m＝∫∫NNTdx1dx2 (16)
N表示全局插值函数矩阵；
得到所有的单元系数矩阵后，根据节点编号和单元编号对应的方式构造全局系数矩阵M和K，从而全局运动控制方程为：
对x2正向和负向扩展到无穷远的区域使用无限元方法进行分析，几何插值函数仍使用双线性插值函数，场插值函数使用衰减型插值函数进行插值离散；x2正向场插值函数为：
公式(19)～(22)中，-1≤ξ≤1,0≤η≤∞，ηi和L为决定衰减程度的任意常数，可以改变ηi/L的值来设定不同的衰减程度；
x2负方向场插值函数为：
公式(23)～(26)中，-1≤ξ≤1,-∞≤η≤0；
无限元方法中场插值函数不同与无限元方法中的插值函数，因此权重函数和概率密度函数的离散放表达式也不同，x2正向表示为：
N′r(x)、N′s(x)为x2正向权函数w和概率密度函数p在某个单元对应的衰减型插值函数；
x2负向表示为：
N″r(x)、N″s(x)为x2负向权函数w和概率密度函数p在某个单元对应的衰减型插值函数；
无限域中的单元系数矩阵的表达式与公式(16)、(17)相似，只需把有限元法中的全局插值函数矩阵N换为无限元方法中的全局插值矩阵N′和N″；将有限域和无限域耦合，构造全局系数矩阵时同样按照节点编号和单元编号一致的原则，按照从左到右，从上到下的原则进行排序；为了构造时逻辑清晰，所以先对正向无限域进行组装，然后是有限元域，最后是负向无限域，由上到下获得全局系数矩阵M和K。
4.如权利要求1所述的Fokker-Planck方程的无限元求解方法，其特征在于，步骤3具体包括：
借助无限元方法，得到每个部分的单元系数矩阵，然后结合算子分裂方法，将单元系数k分裂为k1和k2两项；公式(15)可以写成
其中k1与对流项相关，k2与扩散项相关。k1和k2具体的表达式如下
根据同样的矩阵构原理，得到全局系数矩阵K1和K2，全局运动控制方程可以写为：
。
5.如权利要求1所述的Fokker-Planck方程的无限元求解方法，其特征在于，步骤4具体包括：
使用算子分裂方法求解时，在时间区域[0,T]内，将其划分为n个Δt时间步长；在第i个时间步长内，使用算子分裂连续求解如下两个方程：
在求解式(35)时，以q0＝pi作为初值计算一个时间步长Δt，求得一个解qi+1；在求式(36)时，以p0＝qi+1作为初值计算一个时间步长Δt，求得ti+1时刻的概率密度值pi+1；按照以上求解步骤循环求解，最终可以求得每个时刻的概率密度函数值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311289338.9A CN117350112A (zh) | 2023-10-08 | 2023-10-08 | Fokker-Planck方程的无限元求解方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311289338.9A CN117350112A (zh) | 2023-10-08 | 2023-10-08 | Fokker-Planck方程的无限元求解方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117350112A true CN117350112A (zh) | 2024-01-05 |
Family
ID=89360582
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311289338.9A Pending CN117350112A (zh) | 2023-10-08 | 2023-10-08 | Fokker-Planck方程的无限元求解方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117350112A (zh) |
-
2023
- 2023-10-08 CN CN202311289338.9A patent/CN117350112A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wedig | Dynamic stability of beams under axial forces–Lyapunov exponents for general fluctuating loads | |
| CN108875195B (zh) | 一种考虑接触的三维力学随机振动仿真模拟方法 | |
| Boely et al. | Identification of a non-linear F/A-18 model by the use of fuzzy logic and neural network methods | |
| Vytla et al. | Multi-objective aerodynamic shape optimization of high speed train nose using adaptive surrogate model | |
| CN104504189B (zh) | 随机激励下大规模结构设计方法 | |
| Lichota et al. | D-optimal simultaneous multistep excitations for aircraft parameter estimation | |
| CN104101344A (zh) | 基于粒子群小波网络的mems陀螺随机误差补偿方法 | |
| Wintzer et al. | Optimization and adjoint-based CFD for the conceptual design of low sonic boom aircraft | |
| Rong et al. | Transfer matrix method for dynamics modeling and independent modal space vibration control design of linear hybrid multibody system | |
| Small et al. | A UAV case study with set‐based design | |
| CN115618498A (zh) | 一种飞行器跨流域流场的预测方法、装置、设备及介质 | |
| Mayes et al. | Predicting system response at unmeasured locations | |
| CN117350112A (zh) | Fokker-Planck方程的无限元求解方法 | |
| Nagawkar et al. | Applications of polynomial chaos-based cokriging to aerodynamic design optimization benchmark problems | |
| Dang et al. | Banded null basis and ADMM for embedded MPC | |
| Cooper et al. | Non-intrusive polynomial chaos for efficient uncertainty analysis in parametric roll simulations | |
| Guan et al. | Numerical simulations of sloshing and suppressing sloshing using the optimization technology method | |
| Chin et al. | Analytical modelling of structure-borne sound transmission through I-junction using Chebyshev-Ritz method on cascaded rectangular plate–cavity system | |
| Abraham | Sonic boom loudness reduction through localized supersonic aircraft equivalent-area changes | |
| CN112001100B (zh) | 一种三维地震波场se-ibe耦合模拟方法 | |
| CN108072899B (zh) | 间断Galerkin有限元地震数值模拟算法的自适应实现方法 | |
| Kalashnik | Generation of internal gravity waves by vortex disturbances in a shear flow | |
| Holst | Numerical computation of transonic flow governed by the full-potential equation | |
| CN103984871A (zh) | 一种高精度非线性系统状态估计方法 | |
| Amargianitakis et al. | Generation of noise exposure contours for eVTOL aircraft including transition |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |