王文杰 封建湖

(長安大學(xué)理學(xué)院數(shù)學(xué)系, 西安 710064)

引言

FPK方程于20世紀(jì)初由Fokker與Planck首先提出,并應(yīng)用于研究量子物理問題.20世紀(jì)30年代初Kolmogorov將它一般化與抽象化.不久,Andronov等將它應(yīng)用于研究一般動態(tài)系統(tǒng).50年代Stratonovich將它應(yīng)用于研究電子工程問題.50年代末Chuang與Kazda將它應(yīng)用于研究非線性控制系統(tǒng).60年代初Ariaratnam、Lyon、Smith、Caughey與Dienes、Crandall等將它應(yīng)用于研究非線性隨機(jī)振動問題.

利用轉(zhuǎn)移概率密度求解系統(tǒng)的各種響應(yīng)統(tǒng)計(jì)量,可對該系統(tǒng)的響應(yīng)和可靠性等做定性分析.但是,目前只對一些特殊的非線性隨機(jī)動力系統(tǒng),才能得到其FPK方程的精確解[1-3].基于此,許多學(xué)者致力于FPK方程數(shù)值解法的研究,其中有代表性的主要有：有限元法,有限差分法, 路徑積分法,等價(jià)線性化方法,高斯閉合法,攝動法,Gram-Charlier展開法,等價(jià)非線性系統(tǒng)法,隨機(jī)平均法等.

但這些方法各有其局限性,例如：有限元法一般都很大,而且計(jì)算得到的尾部概率密度不太準(zhǔn)確.等價(jià)線性化方法和高斯閉合法對強(qiáng)非線性[4-8]或具有隨機(jī)參激的系統(tǒng)就不適用,因?yàn)榇藭r(shí)系統(tǒng)響應(yīng)的穩(wěn)態(tài)概率密度往往是非高斯型的.攝動法只適用于弱非線性系統(tǒng).Gram-Charlier展開法可能導(dǎo)致概率密度為負(fù)值的情形.等價(jià)非線性系統(tǒng)法要求兩個(gè)非線性系統(tǒng)性質(zhì)很接近.隨機(jī)平均法只適用于弱阻尼和弱激勵的情形.

而加權(quán)本質(zhì)無振蕩WENO方法,是近年來廣泛流行的一種高分辨率數(shù)值方法,用于解決以對流為主的對流擴(kuò)散方程,特別是雙曲守恒律方程.WENO方法是在ENO方法的基礎(chǔ)上采用加權(quán)思想構(gòu)造的,用于求解包含激波、稀疏博以及接觸間斷等復(fù)雜結(jié)構(gòu)的流體問題.WENO方法性能更穩(wěn)定,對定常問題收斂性更好,它能夠保證在解的光滑區(qū)精度更高,在解的間斷去保持陡峭的間斷過度和本質(zhì)無振蕩性質(zhì).WENO格式最初是在1994年由Liu,Osher和Chan提出,不同于ENO格式單一的模板選取,使用所有候選模板的凸組合,構(gòu)造了一個(gè)三階有限體積WENO格式,系統(tǒng)地討論了WENO方法的構(gòu)造過程和理論分析.隨后在1996年Jiang和Shu在多維空間上給出了一個(gè)可以構(gòu)造任意精度的有限差分WENO[9]格式的框架,并在其中設(shè)計(jì)了一個(gè)五階精度的WENO格式,提出了光滑因子和非線性權(quán)構(gòu)造的基本框架.至今,五階精度的WENO格式是使用最廣泛的,用來求解雙曲守恒律方程WENO重構(gòu)過程的標(biāo)準(zhǔn)格式.

本文針對上述問題,以及結(jié)合三階TVD Runge-Kutta方法,提出了FPK方程的一種新的TVD Runge-Kutta Weno型差分方法,可以較好地解決其他方法所具有的不足,得到比較準(zhǔn)確的概率密度函數(shù).

1 具有隨機(jī)外激非線性動力系統(tǒng)地FPK方程

對于非線性隨機(jī)動力系統(tǒng),其對應(yīng)的FPK方程一般具有如下的形式：

(1)

FPK方程(1)是一個(gè)拋物型變系數(shù)偏微分方程,其描述了擴(kuò)散過程的轉(zhuǎn)移概率密度p(x,t|x0,t0)的進(jìn)化或流動,其中ai,bij為對應(yīng)的漂移和擴(kuò)散系數(shù).

當(dāng)ai,bij均不顯含時(shí)間t時(shí),則FPK方程(1)可轉(zhuǎn)化為：

(2)

FPK方程(2)是一個(gè)橢圓型變系數(shù)偏微分方程,常稱為簡化或平穩(wěn)FPK方程.其解將是平穩(wěn)概率密度p(x).要唯一確定FPK方程(1)的解,還需要初始條件和邊界條件.

在本文中,采用初始條件：

p(x,t|x0,t0)=δ(x-x0),t=t0

(3)

其表示在t=t0時(shí)刻,系統(tǒng)以概率1處于初始狀態(tài)x0.

無窮邊界條件：

(4)

在隨機(jī)振動理論中,FPK方程(1)與(2)常用來預(yù)測非線性隨機(jī)動力系統(tǒng)的響應(yīng).因而,FPK方程解的精確程度對于可靠性分析起著至關(guān)重要的作用.

2 FPK方程的有限差分法

2.1 TVD Runge-Kutta格式

(5)

其中L(u)為f(u)x的逼近.對于半離散格式方程(5),標(biāo)準(zhǔn)Runge-Kutta時(shí)間離散是基于線性穩(wěn)定性條件來實(shí)現(xiàn)格式的穩(wěn)定性,此時(shí)CFL可以取較大的值.但是對于非線性方程,CFL必須很小才能保證格式的穩(wěn)定性,由于高階空間離散和低階Runge-Kutta時(shí)間離散相互耦合,CFL必須大大低于線性穩(wěn)定性要求,才會保證整個(gè)離散格式具有高階精度和無振蕩.本文中采用了具有TVD性質(zhì)的三階Runge-Kutta時(shí)間離散格式[7]：

u(1)=un+ΔtL(un)

其中un是第n時(shí)間層的守恒量,L(u) 是空間離散后的差分算子.

2.2 五階WENO格式

設(shè)f(x)為一函數(shù),則f′(xi)可以利用五階WENO格式表示為：

(6)

其中:

IS0=13(a-b)2+3(a-3b)2,

IS1=13(b-c)2+3(b+c)2,

IS2=13(c-d)2+3(3c-d)2,

Δ+fk=fk+1-fk,Δ-fk=fk-fk-1

其中,ε為小量.

同樣,

而對于f″(xi)利用四階中心差分格式：

(7)

其中fi表示函數(shù)f(x)在點(diǎn)xi處的值,Δxi=xi+1-xi.

2.3 FPK方程的TVD Runge-Kutta WENO型差分格式

對于一維FPK方程(1),將微分的三階TVD Runge-Kutta格式和微分方程的五階WENO格式相結(jié)合,即可得到一維FPK方程的TVD Runge-Kutta WENO型差分格式.

一維問題的TVD Runge-Kutta WENO型差分格式：

(8)

其中：

IS0=13(a-b)2+3(a-3b)2,

IS1=13(b-c)2+3(b+c)2,

IS2=13(c-d)2+3(3c-d)2,

其中,ε為小量.

3 數(shù)值算例

算例1考慮如下受隨機(jī)外激的單自由度非線性系統(tǒng)：

(9)

其中ε=0.1 為常數(shù),標(biāo)識系統(tǒng)的非線性強(qiáng)度,W(t)是一零均值高斯白噪聲,其相關(guān)系數(shù)為：E(W(t)W(t+τ))=2πS0δ(τ) ,S0為W(t)的譜密度,δ(τ)為Dirac函數(shù).

(10)

容易解得：

(11)

其中C為歸一化常數(shù).其精確解與本文方法的數(shù)值解對比如圖1.

圖1 精確解與本文數(shù)值解Fig.1 Exact solution and numerical solution in this paper

算例2考慮如下高斯白噪聲參激和外激聯(lián)合作用下的非線性振子：

X″+2αX′(1+γ1W1(t))+βX′(X2+X′2/ω2)+

ω2X(1+γ2W2(t))=γ3W3(t)

(12)

其中α,β為常數(shù),ω為正常數(shù),Wi(t)(i=1,2,3)為零均值的物理高斯白噪聲,且相互獨(dú)立滿足E[Wi(t)Wi(t+τ)]=δi(τ)(i=1,2,3),δ為Dirac函數(shù),γi(i=1,2,3) 表示各噪聲強(qiáng)度.

圖2 位移y1和速度y2 的邊緣概率密度Fig.2 Marginal probability density of the displacement y1 and the velocity y2

4 結(jié)論

本文將三階TVD Runge-Kutta方法和五階WENO格式相結(jié)合,得到TVD Runge-Kutta WENO型差分格式,并將其成功地應(yīng)用到隨機(jī)外激作用下的非線性動力系統(tǒng),獲得了FPK方程的有限差分?jǐn)?shù)值解,并與其解析解進(jìn)行了比較.表明該方法的有效性以及可行性,并克服一般有限差分法的缺點(diǎn),可以準(zhǔn)確地獲得更小的尾部概率密度,且無振蕩,這對于可靠性分析至關(guān)重要.

1 Lin Y K, Cai G Q. Probabilistic structural dynamics: advanced theory and applications. New York: MeGraw-Hill, 1995

2 朱位秋. 非線性隨機(jī)動力學(xué)與控制——Hamilton理論體系框架. 北京:科學(xué)出版社, 2003 (Zhu W Q. Nonlinear stochastic dynamics and controls——frame of Hamilton theory. Beijing: Science Press, 2003 (in Chinese))

3 趙超櫻,潭維翰,郭志奇. 由菲簡并光學(xué)參量放大系統(tǒng)獲得壓縮態(tài)光所滿足的Fokker-Planck方程及其解. 物理學(xué)報(bào), 2003,52:2694 (Zhao C Y, Tan W H, Guo Z Q. The solution of the Fokker-Planck equation of non-degenerate parametric amplification system for generation of squeezed light.ActaPhysicaSinica, 2003,52:2694 (in Chinese))

4 王平,楊新娥,宋小會. 具有含時(shí)平方反比項(xiàng)的諧振子的路徑積分求解. 物理學(xué)報(bào),2003,52:2957 (Wang P, Yang X E, Song X H. Exact solution for a harmonic oscillator with a time-dependent inverse square potential by path-integral.ActaPhysicaSinica, 2003,52:2957 (in Chinese))

5 徐偉,賀群,戎海武等. Duffing-Van der Pol振子隨機(jī)分叉的全局分析. 物理學(xué)報(bào),2003,52:1365 (Xu W, He Q, Rong H W,et al. Global analysis of stochastic bifurcation in a Duffing-van der Pol system.ActaPhysicaSinica, 2003,52:1365 (in Chinese))

6 孫中奎,徐偉,楊曉麗. 求解強(qiáng)非線性動力系統(tǒng)響應(yīng)的一種新方法. 動力學(xué)與控制學(xué)報(bào),2005(2):29～35 (Sun Z K, Xu W, Yang X L. A new analytic approximate technique for strongly nonlinear dynamic systems.JournalofDynamicsandControl, 2005,2(3):29～35 (in Chinese))

7 Wojtkiewicz S F, Bergman L A, Asce M. Numerical solution of high dimensional Fokker-Planck equations. 8th ASCE Specialty Conference on Probabilistic Mechanics and Structural Reliability, 2000

8 Sun Y, Kumar M. Numerical solution of high dimensional stationary Fokker-Planck equations via tensor decomposition and Chebyshev spectral differentiation.Computers&MathematicswithApplications, 2014,67(10):1960～1977

9 劉儒勛,舒其望. 計(jì)算流體力學(xué)的若干新方法. 北京:科學(xué)出版社,2003:42～106 (Liu R X,Shu Q W. Some new methods for computational fluid dynamics. Beijing:Science Press, 2003:42～106 ( in Chinese))

10 May 2016,revised 19 August 2017.