复杂科学与工程问题仿真的隐式微积分建模

摘要： 针对现代科学与工程仿真遇到愈来愈多难以用经典微积分建模方法描述的复杂问题，在理论研究和工程实践中提出各种含有多个经验参数的唯象偏微分方程模型，或直接采用统计模型来描述.这些模型的物理意义不是很清楚且参数多，其中部分人为参数缺乏物理意义.因此，利用描述问题的基本解或统计分布构造隐式微积分控制方程.这里“隐式”是指可以不需要或难以推导出该控制方程的显式微积分表达式.该方法仅需微积分控制方程的基本解和相应的边界条件就可以进行数值仿真计算.称该方法为隐式微积分方程建模.考虑多相软物质热传导的幂律行为，采用分数阶里斯（Riesz）势核函数为基本解构造稳态问题的隐式分数阶微积分方程模型并进行数值验证.此外，以列维（Lévy）稳态统计分布的概率密度函数为基本解，构造反常扩散现象的隐式分数阶微积分方程模型.该研究的主要数值计算技术基于径向基函数的配点方法.

关键词： 隐式微积分方程建模； 唯象模型； 统计模型； 基本解； 经验参数

中图分类号： O39；O241.8文献标志码： AAbstract： As to a growing number of complex scientific and engineering problems which are not easy to be described by classical calculus modeling methodology， a variety of phenomenological partial differential equation models including multiple empirical parameters have been proposed in theoretical research and engineering practice. In some cases， the statistical models are even used to substitute for the calculus models. These models are not clearly interpreted in physics and require more parameters in which the artificial parameters have no physical significance. Therefore， the fundamental solution or statistical distribution which can describe the problem is employed to construct the implicit calculus governing equation. It is noted that “implicit” in the study suggests that the explicit calculus expression of this governing equation is not required or difficult to derive. The fundamental solution of calculus governing equation and corresponding boundary conditions are sufficient to perform numerical simulation. This strategy is called the implicit calculus equation modeling. Considering the power law behaviors of heat conduction in multiple phase soft materials， the kernel function of fractional Riesz potential is used as the fundamental solution to build the implicit fractional calculus equation model for steadystate problems. The numerical experiments verify the model. In addition， the statistical density function of Lévy stable statistical distribution is used as the fundamental solution to build the implicit calculus equation of fractional order to describe anomalous diffusion. The major numerical technique in the research is the radial basis function based collocation methods.

Key words： implicit calculus equation modeling； phenomenological model； statistical model； fundamental solution； empirical parameter

0引言

微积分是现代数学和古典数学的分水岭，数学的发展和应用自此发生根本性变化.[1]经典的微积分方程建模方法在力学、声学、电磁学、热传输和扩散理论中，甚至在现代量子力学和相对论中取得巨大成功.然而，社会学家、经济学家、物理学家和力学家也发现愈来愈多难以用经典微积分方程建模的所谓“反常”现象[23]，如在扩散和耗散中广泛观察到的幂律现象[34]以及非高斯非马尔科夫过程[56]等.

非线性微分方程模型是描述复杂物理过程的常用方法，已得到充分研究，其基本思路是假设线性力学本构关系或物理定律中的系数是依赖应变变量的.目前，复杂问题的非线性模型愈加复杂，参数很多.例如，岩土力学中的热电化力耦合模型需要四十多个参数，这些参数的物理意义和确定本身就是一个很大的问题.[7]

近年来引起广泛关注的分数阶微积分方法是复杂现象建模的另一个有力的数学工具，在一些领域获得引人注目的成功.[2，4]但是，该方法也有其局限性.首先，非常重要的空间分数阶拉普拉斯算子的定义并不统一，有关数值计算也困难重重[2，8]；其次，分数阶导数阶数的物理解释还不成熟.绝大多数分数阶导数模型都是经验模型或唯象模型.[2，4]

由于实际复杂问题的微分方程模型经常难以建立，因此笔者对这些问题就放弃微分方程建模，直接采用统计模型来描述和分析.[6，9]但是，统计模型不能清晰地描述问题的因果性，物理概念和规律经常不很清楚，结果不精细，一些情况下难以满足实际工程的需要.[5，10]

在微分方程数值模拟方面，目前标准做法是先确定控制方程和边界条件，然后采用某种数值方法做仿真计算.相应的反问题则涉及确定边界条件、控制方程参数和边界形状等，但基本上是先有控制微分方程，然后再求数值解.如上所述，建立复杂问题的微分控制方程并不是一个简单的问题.而且，非线性控制方程和分数阶微分控制方程的数值求解也不是一个容易的任务.例如，边界元法利用微分方程的基本解，能够高效高精度地获得数值解.但是，绝大多数非线性模型的基本解很难找到[11]，而现有的分数阶微积分控制方程的基本解又极为复杂，甚至没有显式表达，也不易得到[2].

为解决仿真这些复杂问题的微积分建模难题，本文提出隐式微积分方程建模方法.基本思路是边界元的逆向思维，即不需要知道微积分控制方程的表达式，而是先确定物理问题微积分方程的基本解或通解，相应的微积分方程存在但不一定能够推出其显式表达式.在数值模拟方面，仅需微积分控制方程的基本解和边界条件就可以进行数值仿真计算，得到模型的数值解，不需要从基本解来推导控制方程.这里“隐式”是指控制方程的显式表达式可以不需要或难以推导出来.在具体实施中可以利用描述一类物理问题的广义基本解或统计分布密度函数.

由于基本解和通解一般可表达为径向基函数，因此本文求解隐式微积分方程模型的主要数值技术是基于径向基函数的配点方法.[12]该类方法以距离为基本变量，不依赖于问题的维数，本质上是无网格无数值积分的方法，编程容易，能够计算高维复杂几何形状问题.

本文考察2类应用实例.首先，考虑多相软物质热传导的幂律行为.[23]许多研究表明：分数阶拉普拉斯算子方程可以有效地描述这类幂律行为的物理力学问题[2，4]，但分数阶拉普拉斯算子的数学定义并不统一[2，13]，现有的表达式都很复杂，难以进行数值计算[14].本文以分数阶里斯（Riesz）势的核函数为基本解构造其稳态问题的隐式微积分方程模型，并用基于径向基函数基本解的奇异边界法[12]进行数值验证.第二个实例是用已知的统计密度函数构造隐式微积分方程的基本解.学者和工程师很早以前就注意到很多工程和社会经济问题不能用经典的高斯分布精确描述，而且难以建立相应的微分方程模型.高斯分布只是列维（Lévy）稳态分布的一个特例[2]，近年来的研究发现列维稳态统计分布比高斯分布应用范围大很多，在许多工程问题得到成功的应用[2，1517]，特别是反常扩散行为中快扩散过程的统计建模.本文运用列维密度函数构造反常扩散现象的时间空间隐式微积分方程模型.本文模型比现有模型简单，物理和统计概念清晰.

本文第1节通过多相软物质幂律热传导建模，引进隐式微积分方程建模方法，并采用奇异边界法给出仿真数值结果，然后在第2节给出列维稳态统计分布的隐式微积分方程模型，最后在第3节总结隐式微积分方程建模方法的特点和优势，以及若干有待研究解决的问题.

①证明过程包含在向J Comput Phys投稿的文章“Threedimensional Rieszkernelbased fractional Laplacian equation and its numerical solution”中，作者为陈文和庞国飞1稳态幂律热传导的隐式微积分方程模型分数阶拉普拉斯算子（-Δ）s/2是一种典型的微分积分算子，能够用单参数s（0到2之间任意实数）表征物理力学系统的空间非局部性；作为经典整数阶拉普拉斯算子（s=2）的一般形式，可用于软物质中声波传播的能量耗散[13]、湍流扩散[16]、地下水溶质运移[1819]、分形空间中的电磁场[20]和非局部热传导[2122]等物理力学问题的建模.算子（-Δ）s/2满足傅里叶变换[8]F{（-Δ）s/2u（·）}=‖k‖sF{u（·）}（1）式中：k为频域中的波数.利用傅里叶逆变换直接推导算子的显式表达式很困难，现有的二维和三维分数阶拉普拉斯算子的显式定义不统一.[13，2224]文献中常用的向量积分显式定义与式（1）不符，是一个近似定义，算子的数值离散也较为困难.例如，有限元离散的弱形式含有二重向量积分，具有非局部性，生成的刚度矩阵不再是带状稀疏阵，而是满阵.[14，21]总之，目前尚无统一的且易于数值计算的分数阶拉普拉斯算子定义.

采用隐式微积分建模方法，笔者不考虑分数阶拉普拉斯算子的具体表达形式，而是从其逆算子（分数阶里斯势）出发，直接构造分数阶拉普拉斯算子的基本解.为不失一般性，三维空间中的分数阶里斯势核函数的定义[8]为u*（x，ξ）=1‖x-ξ‖3-s （2）式中：‖x-ξ‖表示点x与ξ之间的欧氏距离；s为分数阶势的阶数.经典整数阶拉普拉斯算子（s=2）的基本解是分数阶的一个特例，u*（x，ξ）=1‖x-ξ‖ （3）以式（2）作为分数阶拉普拉斯算子（-Δ）s/2的基本解.一般物理问题的分数阶拉普拉斯的阶数s是从1到2之间的实数.可以证明，这样定义的分数阶拉普拉斯算子满足傅里叶变换定义.①

复杂介质往往存在不连续性，如裂纹和孔洞，导致不连续点上的偏导数失去物理意义.经典整数阶导数的微积分方程模型不再适用于描述这类复杂介质中的热传导.[2122]分数阶拉普拉斯方程能够较精确地描述这类幂律（非局部）热传导行为，其稳态方程为-（-Δ）s/2u（x）=0，s∈（1，2]，x∈ΩR3 （4）式中：u为无量纲化的温度函数；s表征材料的非局部性，刻画幂律特征；Ω为计算区域，如图1所示的圆柱.圆柱长为6，底面半径为1，圆柱的中心与坐标原点重合.在本项研究中，（-Δ）s/2按式（2）的分数阶里斯势基本解定义，因此就用这个问题验证基本解式（2）定义的分数阶拉普拉斯算子的隐式微积分模型.需要强调的是，这里并不需要知道分数阶拉普拉斯算子的显式表达式.

基于里斯势的分数阶拉普拉斯算子基本解表达式（2），采用奇异边界法[2526]可直接求解稳态方程式（4）和相应的边界条件的稳态热传导问题.奇异边界法是一种边界型径向基函数配点法，以基本解作为插值基函数.该方法假设基本解源点奇异时的源点强度因子存在.本文采用基本解积分平均计算源点强度因子.

为验证奇异边界法，先考察整数阶拉普拉斯方程（s=2）的数值解精度.图2给出精确解和数值解在圆柱中轴上的值.随着边界离散点数的增加，数值解逐渐逼近精确解，可见奇异边界法具有较好的收敛性.

一般情况下并不知道分数阶拉普拉斯方程式（4）的精确解，但可以通过指定与整数阶方程相同的边界条件考察分数阶方程的数值解是否逼近于整数阶方程的精确解（当s趋于2时）.先考察圆柱中轴{（x，y，z）|x=0，y=0，-3≤z≤3}上的温度随式（4）中分数拉普拉斯算子阶数s的变化，数值结果见图3.在完全相同的边界条件下，当s趋于2时，隐式分数阶拉普拉斯方程的解单调趋近于整数阶拉普拉斯方程的解.此外，s越小，材料的非局部性越强，中轴的温度越低.

2基于列维统计分布的非稳态反常扩散问题的隐式微积分方程模型扩散现象广泛存在于自然界和工业界中，是极其重要的物质迁移和输运的物理力学过程.越来越多的研究发现，经典的扩散方程并不能很好地描述湍流，如高温高压下等离子体扩散，金融市场变化，高分子动力学，以及软物质的热传导、扩散和电子输运等反常扩散过程.所谓的反常扩散[19，27]是指不符合菲克（Fick）扩散定律的扩散行为，包含慢扩散（subdiffusion）和快扩散（superdiffusion）2种形式，通常表现出长程的时间空间相关性.近年来的研究发现空间分数阶扩散方程能较好描述反常扩散中的快扩散现象；但时间空间非稳态分数阶方程的显式表达式难以得到或不准确，且难以数值计算.

本节考虑用列维统计分布的密度函数构造非稳态空间分数阶反常扩散方程的基本解，进行隐式微积分方程建模.这不同于第1节所涉及的稳态问题.

以上分析表明：高斯分布是整数阶菲克扩散模型的基本解核函数，一维列维分布是一维问题分数阶快扩散模型基本解的核函数.列维稳态统计分布是经典扩散方程和空间分数阶扩散方程基本解核函数的两类特殊情况.因此，可以用列维稳态统计分布的概率密度函数构造多维分数阶时间空间扩散方程的基本解，并用于建立快扩散过程的隐式微积分建模.由n维s稳态列维分布概率密度函数得到的n维空间分数阶扩散方程基本解为G（x，y，t）=H（t）tn/sL‖x-y‖t1/s （15）这里列维分布是空间分数阶扩散方程基本解的核函数，深刻揭示多维快扩散过程的统计本质和空间相关性.利用隐式微积分方程模型的基本解式（15），可以用试验或观测数据确定扩散过程所对应的列维统计分布中的稳态指标参数s得到基本解，然后根据可测边界上得到的边界条件值进行数值仿真计算，避免显式表达微积分方程模型的很多困难.

3结论

传统的数学物理方程和数值计算方案一般先根据问题的物理特征和理论采用数学微积分方法建立控制方程和边界条件，然后采用数值方法求解这些偏微分或微分积分方程问题.不同于标准的理论建模和数值仿真方案，本文提出的隐式微积分建模思路是先有问题的基本解，然后直接求解问题.微分控制方程表达式本身不再是必需的环节和对象.

隐式微积分建模的基本解或统计分布可以相当广泛，可极大地推广微积分建模的适用范围.例如，不同于传统的先有微分方程模型再寻找基本解的边界元法，可以直接根据问题的物理特征构造不均匀介质的基本解或通解，甚至可以直接构造非线性问题的基本解，而不用考虑微积分方程的表达形式，可将数学力学建模和数值建模更加紧密地结合起来.

此外，隐式微积分建模方法也将微积分建模与统计模型深刻紧密地结合起来，可由复杂问题的统计分布构造确定性的微分方程模型的基本解，建立确定性模型和随机模型内在联系的桥梁.基本解可以理解为物理场中的影响函数或势函数，由此可建立连续介质的隐式微积分建模与微观尺度的分子动力学和介观尺度的耗散粒子动力学的内在联系.

如何根据复杂问题的物理性质或统计分布构造基本解或通解等影响函数仍是有待深入研究的课题.

致谢：本文的第1节和第2节分别得到博士研究生庞国飞和博士傅卓佳的帮助，在此表示感谢.

参考文献：

[1]莫里斯·克莱因. 古今数学思想[M]. 张理京， 译. 上海： 上海科学技术出版社， 2009： 342383.

[2]陈文， 孙洪广， 李西成， 等. 力学与工程问题的分数阶导数建模[M]. 北京： 科学出版社， 2010： 8285.

[3]马红孺， 陆坤权. 软凝聚态物质物理学[J]. 物理， 2000， 29（9）： 561524.

MA Hongru， LU Kunquan. The physics of soft condensed matter[J]. Physics， 2000， 29（9）： 561523.

[4]徐明瑜， 谭文长. 中间过程、临界现象——分数阶算子理论、方法、进展及其在现代力学中的应用[J]. 中国科学： G辑： 物理学力学天文学， 2006， 36（3）： 225238.

XU Mingyu， TAN Wenchang. Intermediate process， critical phenomena： theories， methods， and advances of fractional operator and their applications in modern mechanics[J]. Sci China： Ser G： Phys， Mech & Astron， 2006， 36（3）： 225238.

[5]MANDELBROT B B. The fractal geometry of nature[M]. New York： WH Freeman， 1982： 247272.

[6]METZLER R， KLAFTER J. The random walks guide to anomalous diffusion： a fractional dynamics approach[J]. Phys Rep， 2000（339）： 177.

[7]周创兵， 陈益峰， 姜清辉， 等. 论岩体多场广义耦合及其工程应用[J]. 岩石力学与工程学报， 2008， 28（7） ： 13291340.

ZHOU Chuangbing， CHEN Yifeng， JIANG Qinghui， et al. On generalized multifield coupling for fractured rock masses and its applications to rock engineering[J]. Chin J Rock Mech & Eng， 2008， 28（7）： 13291340.

[8]SAMKO S G， KILBAS A A， MARICHEV O I. Fractional integrals and derivatives： theory and applications[M]. London： Gordon and Breach Science Publishers， 1993： 483532.

[9]SHLESINGER M F. Fractal time and 1/f noise in complex systems[J]. Ann NY Acad Sci， 1987（504）： 214228.

[10]谢和平. 分形岩石力学导论[M]. 北京： 科学出版社， 2005.

[11]KYTHE P K. Fundamental solutions for differential operators and applications[M]. Boston： Birkhauser， 1996： 60136.

[12]CHEN W， FU Z， CHEN C. Recent advances in radial basis function collocation methods[M]. SpringerVerlag， 2013.

[13]CHEN W， HOLM S. Fractional Laplacian timespace models for linear and nonlinear lossy media exhibiting arbitrary frequency powerlaw dependency[J]. J Acoust Soc Am， 2004， 115（4）， 14241430.

[14]ROOP J P. Computational aspects of FEM approximation of fractional advection diffusion equations on boundary domains in R2[J]. J Comput Appl Math， 2006， 193（1）： 243268.

[15]DELCASTILLONEGRETE D， CARRERAS B A， LYNCH V E. Front dynamics in reactiondiffusion systems with Levy flights： a fractional diffusion approach[J]. Phys Rev Lett， 2003： 91（1）： 0183014.

[16]CHEN W. A speculative study of 2/3order fractional Laplacian modeling of turbulence： some thoughts and conjectures[J/OL]. Chaos， 2006（16）. [20140827]. http：//dx.doi.org/10.1063/1.2208452.

[17]CHEN W. Lévy stable distribution and[0，2] power law dependence of acoustic absorption on frequency in various lossy media[J]. Chin Phys Lett， 2005， 22（10）： 26012603.

[18]MEERSCHAERT M M， BENSON D A， BAUMER B. Multidimensional advection and fractional dispersion[J]. Phys Rev E， 1999， 59（5）： 50265028.

[19]PANG G， CHEN W， FU Z. Spacefractional advectiondispersion equations by the Kansa method[EB/OL].（20140720）[20140830]http：//www.sciencedirect.com/science/article/pii/S0021999114005130.

[20]TARASOV V E. Electromagnetic fields on fractals[J]. Mod Phys Lett A， 2006， 21（20）： 15871600.

[21]BOBARU F， DUANPANYA M. The peridynamic formulation for transient heat conduction[J]. Int J Heat & Mass Tran， 2010， 53（19/20）： 40474059.

[22]DELIA M， GUNZBURGER M. The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator[J]. Comput & Math Appl， 2013， 66（7）： 12451260.

[23]TARASOV V E. Fractional vector calculus and fractional Maxwells equations[J]. Ann Phy， 2008， 323（11）： 27562778.

[24]MEERSCHAERTA M M， MORTENSEN J， WHEATCRAFTS W. Fractional vector calculus for fractional advectiondispersion[J]. Physica A： Stat Mech & its Appl， 2006（367）： 181190.

[25]陈文. 奇异边界法： 一个新的、简单、无网格、边界配点数值方法[J]. 固体力学学报， 2009， 30（6）： 592599.

CHEN Wen. Singular boundary method： a novel， simple， meshfree boundary collocation numerical method[J]. Chin J Solid Mech， 2009， 30（6）： 592599.

[26]GU Y， CHEN W， HE X Q. Improved singular boundary method for elasticity problems[J]. Computer and Structures， 2014， 135： 7382.

[27]FU Z， CHEN W， YANG H. Boundary particle method for Laplace transformed time fractional diffusion equations[J]. J Comput Phys， 2013， 235（15）： 5266.

[28]FELLER W. An introduction to probability theory and its applications： Vol 2 [M]. 2nd ed. New York： Wiley， 1971： 574581.

[29]SAICHEV A I， ZASLAVSKY G M. Fractional kinetic equations： solutions and applications[J]. Chaos， 1997， 7（4）： 753764.

[30]LIANG Y， CHEN W. A survey on computing Lévy stable distributions and a new MATLAB toolbox[J].

[19]PANG G， CHEN W， FU Z. Spacefractional advectiondispersion equations by the Kansa method[EB/OL].（20140720）[20140830]http：//www.sciencedirect.com/science/article/pii/S0021999114005130.

[20]TARASOV V E. Electromagnetic fields on fractals[J]. Mod Phys Lett A， 2006， 21（20）： 15871600.

[21]BOBARU F， DUANPANYA M. The peridynamic formulation for transient heat conduction[J]. Int J Heat & Mass Tran， 2010， 53（19/20）： 40474059.

[22]DELIA M， GUNZBURGER M. The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator[J]. Comput & Math Appl， 2013， 66（7）： 12451260.

[23]TARASOV V E. Fractional vector calculus and fractional Maxwells equations[J]. Ann Phy， 2008， 323（11）： 27562778.

[24]MEERSCHAERTA M M， MORTENSEN J， WHEATCRAFTS W. Fractional vector calculus for fractional advectiondispersion[J]. Physica A： Stat Mech & its Appl， 2006（367）： 181190.

[25]陈文. 奇异边界法： 一个新的、简单、无网格、边界配点数值方法[J]. 固体力学学报， 2009， 30（6）： 592599.

CHEN Wen. Singular boundary method： a novel， simple， meshfree boundary collocation numerical method[J]. Chin J Solid Mech， 2009， 30（6）： 592599.

[26]GU Y， CHEN W， HE X Q. Improved singular boundary method for elasticity problems[J]. Computer and Structures， 2014， 135： 7382.

[27]FU Z， CHEN W， YANG H. Boundary particle method for Laplace transformed time fractional diffusion equations[J]. J Comput Phys， 2013， 235（15）： 5266.

[28]FELLER W. An introduction to probability theory and its applications： Vol 2 [M]. 2nd ed. New York： Wiley， 1971： 574581.

[29]SAICHEV A I， ZASLAVSKY G M. Fractional kinetic equations： solutions and applications[J]. Chaos， 1997， 7（4）： 753764.

[30]LIANG Y， CHEN W. A survey on computing Lévy stable distributions and a new MATLAB toolbox[J].

[19]PANG G， CHEN W， FU Z. Spacefractional advectiondispersion equations by the Kansa method[EB/OL].（20140720）[20140830]http：//www.sciencedirect.com/science/article/pii/S0021999114005130.

[20]TARASOV V E. Electromagnetic fields on fractals[J]. Mod Phys Lett A， 2006， 21（20）： 15871600.

[21]BOBARU F， DUANPANYA M. The peridynamic formulation for transient heat conduction[J]. Int J Heat & Mass Tran， 2010， 53（19/20）： 40474059.

[22]DELIA M， GUNZBURGER M. The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator[J]. Comput & Math Appl， 2013， 66（7）： 12451260.

[23]TARASOV V E. Fractional vector calculus and fractional Maxwells equations[J]. Ann Phy， 2008， 323（11）： 27562778.

[24]MEERSCHAERTA M M， MORTENSEN J， WHEATCRAFTS W. Fractional vector calculus for fractional advectiondispersion[J]. Physica A： Stat Mech & its Appl， 2006（367）： 181190.

[25]陈文. 奇异边界法： 一个新的、简单、无网格、边界配点数值方法[J]. 固体力学学报， 2009， 30（6）： 592599.

CHEN Wen. Singular boundary method： a novel， simple， meshfree boundary collocation numerical method[J]. Chin J Solid Mech， 2009， 30（6）： 592599.

[26]GU Y， CHEN W， HE X Q. Improved singular boundary method for elasticity problems[J]. Computer and Structures， 2014， 135： 7382.

[27]FU Z， CHEN W， YANG H. Boundary particle method for Laplace transformed time fractional diffusion equations[J]. J Comput Phys， 2013， 235（15）： 5266.

[28]FELLER W. An introduction to probability theory and its applications： Vol 2 [M]. 2nd ed. New York： Wiley， 1971： 574581.

[29]SAICHEV A I， ZASLAVSKY G M. Fractional kinetic equations： solutions and applications[J]. Chaos， 1997， 7（4）： 753764.

[30]LIANG Y， CHEN W. A survey on computing Lévy stable distributions and a new MATLAB toolbox[J].