復雜科學與工程問題仿真的隱式微積分建模

摘要： 針對現代科學與工程仿真遇到愈來愈多難以用經典微積分建模方法描述的復雜問題，在理論研究和工程實踐中提出各種含有多個經驗參數的唯象偏微分方程模型，或直接采用統計模型來描述.這些模型的物理意義不是很清楚且參數多，其中部分人為參數缺乏物理意義.因此，利用描述問題的基本解或統計分布構造隱式微積分控制方程.這里“隱式”是指可以不需要或難以推導出該控制方程的顯式微積分表達式.該方法僅需微積分控制方程的基本解和相應的邊界條件就可以進行數值仿真計算.稱該方法為隱式微積分方程建模.考慮多相軟物質熱傳導的冪律行為，采用分數階里斯（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].