康紅喜, 張引娣, 蔣 茜
(長(zhǎng)安大學(xué) 理學(xué)院,陜西 西安 710064)
隨著隨機(jī)微分方程在金融、生物、工程技術(shù)等領(lǐng)域應(yīng)用變得越來(lái)越多,對(duì)于隨機(jī)微分方程的數(shù)值解的要求也越來(lái)越高。為了滿足各個(gè)領(lǐng)域的需求,學(xué)者們提出了許多求解隨機(jī)微分方程的數(shù)值方法,文獻(xiàn)[1-8]給出了常見數(shù)值方法相關(guān)的研究成果。
求解微分方程的方法還有Heun方法,但是該方法對(duì)步長(zhǎng)是有要求的,學(xué)者們的研究相對(duì)較少。Heun方法在文獻(xiàn)[9]中被提出來(lái)后,文獻(xiàn)[10-11]研究了該方法收斂性與穩(wěn)定性的條件。另外,求解剛性隨機(jī)微分方程對(duì)步長(zhǎng)的要求更為苛刻,為了解決這一問(wèn)題,文獻(xiàn)[12]提出了一種全隱式的方法,又稱平衡法。文獻(xiàn)[13-15]對(duì)這個(gè)全隱式方法作了進(jìn)一步研究,得到了該方法的漸近穩(wěn)定性、均方穩(wěn)定性以及高階收斂性。本文提出了一個(gè)平衡θ-Heun方法,該方法結(jié)合了平衡方法與文獻(xiàn)[16]中θ-Heun法的優(yōu)點(diǎn)。本文研究該方法的均方穩(wěn)定性和指數(shù)穩(wěn)定性,并比較3種Heun方法的穩(wěn)定性。
考慮d維標(biāo)量自治隨機(jī)微分方程
(1)
其中:t∈[0,T];x∈Rd;函數(shù)f(x)為漂流項(xiàng),函數(shù)g(x)為擴(kuò)散項(xiàng),兩者在[0,T]上都是連續(xù)可測(cè)的,且有E|x0|2<∞;W(t)為一維標(biāo)準(zhǔn)布朗運(yùn)動(dòng)。設(shè)h為步長(zhǎng),當(dāng)t>0,h>0時(shí),其增量ΔW(t)=W(t+h)-W(t)服從正態(tài)分布N(0,h)。本文只研究方程(1)的一個(gè)類型,即帶乘性噪音的隨機(jī)微分方程。另外迭代式xn+1=R(λ,μ,h,ΔWn)xn是本文用數(shù)值方法求解隨機(jī)微分方程得到的。
1.2.1 θ-Heun方法
θ-Heun方法計(jì)算公式[16]為:
xn+1=xn+(1-θ)f(xn)h+θf(wàn)(xn+hf(xn))h+g(xn)ΔWn。
1.2.2 平衡法
平衡法計(jì)算公式為[17]:
xn+1=xn+f(xn)h+g(xn)ΔWn+Cn(xn-xn+1),
其中,Cn=c0(tn,xn)h+c1(tn,xn)|ΔWn|,函數(shù)c0、c1稱為控制函數(shù)。
1.2.3 平衡θ-Heun方法
定義1 令h=T/N,tn=nh,n=0,1,…,N,稱
xn+1=xn+(1-θ)f(xn)h+g(xn)ΔWn+θf(wàn)(xn+hf(xn))h+Cn(xn-xn+1)
(2)
是平衡θ-Heun方法。其中:θ∈[0,1];Cn=c0(tn,xn)h+c1(tn,xn)|ΔWn|,特別地,當(dāng)θ=0時(shí),(2)式就是平衡隱式方法。根據(jù)文獻(xiàn)[15],若要保證系統(tǒng)的穩(wěn)定性,則這里的c0、c1是有界的d×d的實(shí)值矩陣值函數(shù),滿足可逆矩陣
M(t,x)=I+α0c0(t,x)+α1c1(t,x),
|M(t,x)-1|≤K<∞
(3)
E|xn|2≤NE|x0|2e-mnh, ?n≥0
(4)
成立,則稱數(shù)值方法是指數(shù)穩(wěn)定的。
為了能夠更簡(jiǎn)便地研究數(shù)值方法的穩(wěn)定性,引入下面的隨機(jī)微分方程,即
(5)
其中:t∈[0,T];λ,μ∈R;x∈Rd;W(t)為一維標(biāo)準(zhǔn)布朗運(yùn)動(dòng)。
將平衡θ-Heun方法用于(5)式得到如下的差分方程:
xn+1=xn+λhxn+θλ2h2xn+μxnΔWn+Cn(xn-xn+1)
(6)
其中,Cn=c0h+c1|ΔWn|,這里假設(shè)c0、c1為常數(shù),步長(zhǎng)為h=T/N,T=Nt,tn=nh,n=0,1,…,N。
定理1 用平衡θ-Heun方法求解線性隨機(jī)微分方程(5)時(shí),若有:
則方程(5)是MS穩(wěn)定的。
證明首先將差分方程(6)寫成下面的形式:
xn+1=(1+Cn)-1[(1+λh+θλ2h2+μΔWn)+Cn]xn=(1+Cn)-1[(1+λh+θλ2h2+μΔWn)+c0h+c1|ΔWn|]xn
(7)
根據(jù)(3)式有:
(8)
將(7)式、(8)式兩邊同時(shí)取均值可得:
(9)
(10)
其中
2K2E((c1+μ)|ΔWn|)2。
遞推(10)式,要使Sn→0(n→∞),只需
根據(jù)期望的定義及積分性質(zhì)可得:
即
1-2K2(1+λh+θλ2h2+c0h)2。
因此,若要方程(5)是MS穩(wěn)定的,則只需:
下面取特殊的參數(shù),c0=-λ/2,c1=0,當(dāng)方程(5)為一維時(shí),證明平衡θ-Heun方法MS穩(wěn)定的具體條件以及平衡θ-Heun方法的指數(shù)穩(wěn)定性。
q<-(θ2p4+θp3+2θp2+2p),
則平衡θ-Heun方法是MS穩(wěn)定的。其中:p=λh;q=μ2h。
證明由(6)式可得:
(11)
(11)式中取c0=-λ/2,c1=0,令p=λh,q=μ2h,則有:
(12)
對(duì)(12)式兩邊同時(shí)平方取均值,由E(ΔWn)=0,E(ΔWn)2=h可得:
(13)
整理(13)式可得:
q<-(θ2p4+θp3+2θp2+2p)。
此時(shí)平衡θ-Heun方法是MS穩(wěn)定的。
定理3對(duì)于方程(5),當(dāng)c0=-λ/2,c1=0時(shí),平衡θ-Heun方法是指數(shù)穩(wěn)定的,當(dāng)且僅當(dāng)有下面的不等式成立,即
證明將c0=-λ/2,c1=0帶入(6)式,可得:
(14)
對(duì)(14)式兩邊分別平方取均值,由E(ΔWn)=0,E(ΔWn)2=h,可得:
(15)
記
則(15)式經(jīng)過(guò)遞推可得:
(16)
當(dāng)c0=-λ/2,c1=0時(shí),求解方程(5)的平衡θ-Heun方法是指數(shù)穩(wěn)定的,則根據(jù)定義3有:
E|xn|2≤NE|x0|2e- mnh, ?n≥0。
若想得到不等式
則需要找到一個(gè)n,使得:
(17)
又由:
(18)
結(jié)合(16)式、(18)式可得:
E|xn|2≤NE|x0|2e- mnh, ?n≥0,
則當(dāng)c0=-λ/2,c1=0時(shí),平衡θ-Heun方法是指數(shù)穩(wěn)定的。
平衡θ-Heun方法、Heun方法、θ-Heun方法的MS穩(wěn)定性如圖1所示。
圖1 3種方法的MS穩(wěn)定性
從圖1可以看出,平衡θ-Heun方法的穩(wěn)定性最好,且對(duì)給定的3個(gè)步長(zhǎng)都是MS穩(wěn)定的;而Heun方法、θ-Heun方法的穩(wěn)定性較差,只有h=1/4時(shí)才保證了方法的穩(wěn)定性。通過(guò)對(duì)比可知,3種數(shù)值方法中平衡θ-Heun方法對(duì)步長(zhǎng)的要求最低。 下面進(jìn)行理論分析與驗(yàn)證。
(1) 對(duì)于圖1a,設(shè)
θλ3h3+θ2λ4h4+μ2h,
當(dāng)步長(zhǎng)為h=1,h=1/2,h=1/4時(shí),對(duì)應(yīng)的值為:
由定義2可知,對(duì)于3個(gè)步長(zhǎng),平衡θ-Heun方法都滿足MS穩(wěn)定條件。
(2) 對(duì)于圖1b,設(shè)
當(dāng)步長(zhǎng)為h=1,h=1/2,h=1/4時(shí),對(duì)應(yīng)的值為:
由定義2可知,當(dāng)h=1/4時(shí),Heun法才滿足MS穩(wěn)定的條件。
(3) 對(duì)于圖1c,設(shè)
c=(1+λh+θλ2h2)2+μ2h,
當(dāng)步長(zhǎng)為h=1,h=1/2,h=1/4時(shí),對(duì)應(yīng)的值為:
由定義2可知,當(dāng)h=1/4時(shí),θ-Heun方法滿足MS穩(wěn)定的條件。
本文將平衡法的思想與θ-Heun方法結(jié)合起來(lái),得到了平衡θ-Heun方法,研究了該方法用于d維隨機(jī)微分方程(5)系統(tǒng)穩(wěn)定的條件;在給定合適的控制參數(shù)時(shí),對(duì)于一維線性隨機(jī)微分方程,得到了保證平衡θ-Heun方法MS穩(wěn)定的條件以及指數(shù)穩(wěn)定的條件。從最后的數(shù)值算例中可以看出,平衡θ-Heun方法在3種數(shù)值方法中對(duì)于步長(zhǎng)的要求最低。