譚加博,蔡成偉

(北京物資學(xué)院信息學(xué)院,北京 101149)

一類二階微分方程數(shù)值新算法

譚加博,蔡成偉

(北京物資學(xué)院信息學(xué)院,北京 101149)

通過辛Euler數(shù)值算法的構(gòu)造過程的研究,在發(fā)現(xiàn)構(gòu)造辛Euler法格式的同時,還有另外的兩個相似性質(zhì)的算法格式,對這兩種算法格式的研究以及與辛Euler算法的比較,得出數(shù)值算法優(yōu)秀的標(biāo)準(zhǔn)以及優(yōu)秀數(shù)值算法形成過程中需要改進的方向.通過擺方程比較這些數(shù)值算法實際的計算效果,得知辛Euler算法格式優(yōu)于其他兩種相似性質(zhì)的算法格式.

二階微分方程;時間代價;辛Euler法;新數(shù)值算法;擺方程

0 引言

通過對一階微分方程的Euler算法、改進Euler算法、隱Euler算法、Gauss中點算法以及對二階微分方程的數(shù)值解法辛Euler法的研究,發(fā)現(xiàn)在一階微分方程這組算法中,用于計算二階微分方程的數(shù)值算法的辛Euler算法是由Euler算法和隱Euler算法組合而成.于是在研究辛Euler法的同一理論平臺上,筆者提出了一個由Euler算法和Gauss中點算法組合或者是由Euler算法和改進Euler算法組合而成的新算法,并對這兩種新組合方法在求解二階微分方程中運用,與目前流行的辛Euler法做比較,從而發(fā)現(xiàn)所給新算法的優(yōu)缺點.

1 4種基本數(shù)值算法的格式推導(dǎo)

其中u0是給定的初值,這就是一階常微分方程的初值問題.為使問題(1)的解存在、唯一且連續(xù)依賴初值即初值u0問題(1)適定,還必須對右端f(t,u)加適當(dāng)限制,通常要求f關(guān)于u滿足Lipschitz條件:即存在常數(shù)L,使

對所有t∈[0,T]和u1,u2∈(-∞,∞)成立.在常微分方程數(shù)值解法中,對于常微分方程初值問題的數(shù)值解法都要滿足上述的Lipschitz條件[1].

1.1 Euler法

最簡單的數(shù)值解法是Euler法.將區(qū)間[0,T]作N等分,小區(qū)間的長度h=T/N稱為步長,點列tn=nh(n=0,1,2,…,N)稱為節(jié)點,t0=0由已知初值u(t0)=u0,可算出u(t)在t=t0的導(dǎo)數(shù)值u′(t0)=f(t0,u(t0))=f(t0,u0).利用Taylor展式

其中ζ∈(t0,t1),并略去二階小量R0,得u1=u0+hf(t0,u0).u1就是u(t1)的近似值.利用u1又可算出u2,如此下去可算出u在所有節(jié)點上的值,一般遞推公式為

這就是Euler法[2].

1.2 改進Euler法

下面對Euler法運用數(shù)值積分法推導(dǎo).將問題(1)寫成等價的積分形式

特別地

用左矩形公式近似右端積分,并用u1替代u(t1),即得u1=u0+hf(t0,u0),這就是Euler法(式(5)).我們也可用梯形公式近似上述積分,仍用u1替代u(t1),得

稱之為改進Euler法[2].

1.3 隱式Euler法

隱式Euler法[3]同Euler法的數(shù)值解法相似,將區(qū)間[0,T]作N等分,小區(qū)間的長度h=T/N稱為步長,點列tn=nh(n=0, 1,2,…,N)稱為節(jié)點,t0=0是已知初值u(t0)=u0.在Euler法中,數(shù)值計算的迭代格式中運用初始值點(t0,u0)進行格式迭代計算.在改進Euler法中,數(shù)值計算的迭代格式中運用初始值點(t0,u0)和所求值點(t1,u1)進行格式迭代計算.根據(jù)上面的Euler法和改進Euler法的格式推導(dǎo)過程可以進行相似的演變.將Euler法中迭代點(t0,u0)用所求值點(t1,u1)代替,因此有數(shù)值計算的迭代格式

1.4 Gauss中點法

對于問題(1)上面3種數(shù)值算法都可求解,但是,這樣的求解結(jié)果并沒有得到微分方程數(shù)值解研究者的滿足.我們需要得到更好的數(shù)值算法.Gauss在對上述的3種數(shù)值算法反復(fù)研究后,再次改造Euler算法的算式格式.將運用迭代格式中相鄰2點的中點作為算式的迭代值,即取t=0初始點(t0,u0)與后一點(t1,u1),將這2點的中點((t0+t1)/2,(u0+u1)/2)作為Euler法格式中的迭代點.即,這個中點迭代格式為

在一般情況下,將區(qū)間[0,T]作N等分,小區(qū)間的長度h=T/N稱為步長,點列tn=nh(n=0,1,2,…,N)稱為節(jié)點,迭代格式為

這個迭代格式就稱為Gauss中點法[3].

2 辛Euler算法推導(dǎo)以及誤差估計、局部穩(wěn)定性

通過4個基本數(shù)值算法,得到一階微分方程的數(shù)值求解.但是,對于二階微分方程

的數(shù)值解,通過一階微分方程的數(shù)值算法思想知道可將二階微分方程化為兩個一階微分方程來求數(shù)值解.在此假設(shè)令y′= w(x,y),對二階微分方程(10)可以化為

通過求解一階微分方程組來求解這個二階微分方程.應(yīng)用上面的4種基本數(shù)值算法來求解一階微分方程組.

通過上述的微分方程轉(zhuǎn)化后,將Euler法和隱Euler法應(yīng)用到二階微分方程化成的兩個一階微分方程中.將微分方程組(11)中的第一個微分方程用隱Euler法求數(shù)值解,第二個微分方程用Euler法求數(shù)值解,然后將兩個數(shù)值解格式合并在一起后的數(shù)值算法格式,稱為辛Euler法[1].由于辛Euler算法是在局部穩(wěn)定性算法的基礎(chǔ)上組合而成,同時微分方程也是滿足Lipschitz條件,所以辛Euler數(shù)值算法也是局部穩(wěn)定.對于上面的二階微分方程,假設(shè)在數(shù)值算法中,將Euler算法和隱Euler算法的迭代步長相等為h代入,得數(shù)值算法格式推導(dǎo)過程如下.

格式(12)為辛Euler法的通式.

將二階微分方程的真解,應(yīng)用泰勒展開式展開后有

可以看出,辛Euler法的局部截斷誤差en的階次為O(h3).辛Euler法的整體截斷誤差εn的階次為O(h2).

通過辛Euler算法的誤差估計可知,辛Euler法也是收斂的.在此延續(xù)前面的假設(shè)條件,對于辛Euler法的格式(12)的局部穩(wěn)定性研究如下

在這里L(fēng),h,e0都是常數(shù),所以可以知道辛Euler法是局部穩(wěn)定的.同時,辛Euler法是顯示格式.

3 另外兩種數(shù)值算法以及誤差估計、局部穩(wěn)定性

在辛Euler法的推導(dǎo)中,可以知道,在辛Euler法產(chǎn)生的算法理論平臺上,除了有Euler法與隱Euler法可以組成顯示迭代格式外,還有另外的兩種基本算法的組合方式也能形成顯示迭代格式.在這兩種組合方式中,一個算法是由Euler法與改進Euler法組合而成;一個算法是由Euler法與Gauss中點法組合而成.這3種數(shù)值算法都是同一個算法理論平臺中的顯示格式.從前面的辛Euler法中可以知道,這兩種算法也是收斂的.

3.1 第一種數(shù)值算法

在二階微分方程數(shù)值求解中,將Euler法與改進Euler法組合而成.對于二階微分方程(10)問題的求解,將其轉(zhuǎn)化為一階微分方程組(11)的形式.在此借用辛Euler法中的假設(shè)條件,將第一個一階微分方程用改進Euler法格式,第二個一階微分方程用Euler法格式,故有此數(shù)值算法的格式推導(dǎo)如下

這個格式就是第一種數(shù)值算法的格式.第一種數(shù)值算法格式與辛Euler法的格式在最后一部分上系數(shù)不同.

將二階微分方程的真解,應(yīng)用泰勒展開式展開后有

可以看出,第一種數(shù)值算法的局部截斷誤差en的階次為O(h3).通過改進Euler法的結(jié)論可以知道第一種數(shù)值算法的整體截斷誤差εn的階次為O(h2).

通過上面的第一種數(shù)值算法的格式與二階微分方程真解的展開格式相比可見,前者更加接近真解.第一種數(shù)值算法的格式與真解的展開式的前面3部分值完全相同.在它們誤差相當(dāng)?shù)那闆r下,從數(shù)值算法的格式中可見,第一種數(shù)值算法在構(gòu)造的迭代格式上比辛Euler法要更精確些.

第一種數(shù)值算法在此延續(xù)前面的假設(shè)條件,對于第一種數(shù)值算法的格式(13)的局部穩(wěn)定性[3]研究如下

在這里L(fēng),h,e0都是常數(shù),所以第一種數(shù)值算法格式是局部穩(wěn)定的.

3.2 第二種數(shù)值算法

在二階微分方程數(shù)值求解中,將Euler法與Gauss中點法組合而成.對于二階微分方程(10)問題的求解,將其轉(zhuǎn)化為一階微分方程組(11)的形式.在此借用辛Euler法中的假設(shè)條件,將第一個一階微分方程用Euler法格式,第二個一階微分方程用Gauss中點法格式,有這個數(shù)值算法的格式推導(dǎo)如下

將真解的展開式同第二種數(shù)值算法的迭代格式的轉(zhuǎn)化形式(18)相比較,第二種數(shù)值算法格式的局部截斷誤差en的階次為O(h3),整體截斷誤差εn的階次為O(h2).

第二種數(shù)值算法的迭代格式在形式上比辛Euler法的迭代格式更加復(fù)雜,所以在數(shù)值計算時需要的初始條件更多.與辛Euler法相比,在誤差上相當(dāng).從數(shù)值算法的構(gòu)造迭代格式上來看,辛Euler法要好于第二種數(shù)值算法.

第二種數(shù)值算法在此延續(xù)前面的假設(shè)條件,對于第二種數(shù)值算法的格式(14)的局部穩(wěn)定性研究如下

這里L(fēng),h,e0都是常數(shù).從上面的不等式可知,第二種數(shù)值算法格式的局部截斷誤差en是有上限的,且由前面的兩步迭代截斷誤差決定,所以第二種數(shù)值算法格式具有局部穩(wěn)定性[3].

4 3種二階微分方程數(shù)值算法的實驗比較

衡量算法的時間代價也是評價一種算法的指標(biāo),所以,對于數(shù)值算法的研究,必須通過計算機的實際檢驗.同時,對于算法的整體穩(wěn)定性也是通過算法的實際計算得到的結(jié)果.

對于二階微分方程的數(shù)值算法中選擇的辛Euler算法,以及在同一個算法的理論平臺中,還有另外兩種數(shù)值算法,它們都是顯示算法格式.同時,通過算法的實驗,也能夠很好地直觀地感覺算法的精確性、穩(wěn)定性.為了算法實驗獲得更直觀的感受,將算法通過擺方程來實驗,擺方程在實際的計算中有很好的視覺效果.在算法實驗前,先假設(shè)擺方程的初始問題如下

在此,假設(shè)q′=p,則擺方程化成一階微分方程組有

在擺方程的兩個一階微分方程的數(shù)值計算中,假設(shè)兩個微分方程的數(shù)值計算步長相同,設(shè)為h.迭代計算的總步長數(shù),設(shè)為N,則n=0,1,…,N.其中,將擺方程問題的初始值是t=0的值記為p0=0,q0=0.5.在算法實驗中,為了知道算法的運算時間代價,將算法的運算完成的時間記為T.

4.1 辛Euler法擺方程實驗算式格式

對于擺方程初始值問題(19)在辛Euler法中的實驗,由辛Euler法的格式和初始值的擺方程問題有實驗的算法迭代格式

4.2 第一種數(shù)值算法擺方程算式格式

對于擺方程初始值問題(19)在第一種數(shù)值算法中的實驗,由第一種數(shù)值算法的格式和初始值的擺方程問題有實驗的算法迭代格式

4.3 第二種數(shù)值算法擺方程算式格式

對于擺方程初始值問題(19)在第二種數(shù)值算法中的實驗,由第二種數(shù)值算法的格式和初始值的擺方程問題有實驗的算法迭代格式

下面對3個算法的擺方程問題數(shù)值計算運用不同的步長h,不同的迭代步長數(shù)N編程,并將運行結(jié)果以圖像的形式顯示[4].

4.4 3種數(shù)值算法整體穩(wěn)定性比較分析

在h=0.1,N=400時3種數(shù)值算法的計算結(jié)果如圖1.可見,第二種數(shù)值算法的整體穩(wěn)定性是不好的,在迭代不到400步就向負(fù)方向飛出.在h=0.1,N=900時3種數(shù)值算法的計算結(jié)果如圖2,可知第一種數(shù)值算法在迭代接近900步時,算法的迭代數(shù)值就飛離了,所以第一種數(shù)值算法的整體穩(wěn)定性也不是很好.從3種數(shù)值算法的擺方程的7個計算結(jié)果圖中能夠很明顯地看出,第一種數(shù)值算法和第二種數(shù)值算法的迭代運算值是向外擴散的,而辛Euler法的迭代運算值一直很穩(wěn)定地保持在一條線帶上,所以辛Euler法的整體穩(wěn)定性比另外兩種數(shù)值算法好.

圖7 3種算法計算擺方程 h=0.001,N=100 000Fig.7 Three methods of pendulum equation h=0.001,N=100 000

4.5 3種數(shù)值算法的時間代價

在實驗中,通過在同一臺計算機上進行3種數(shù)值算法計算同一個擺方程的實驗,同時在相同的計算機環(huán)境下平衡實驗的系統(tǒng)誤差,去除實驗的偶然誤差,對算法運算的計時知道,在步數(shù)N<10 000時,都能夠很快地完成,對于步數(shù)N>10 000時,都需要一段時間才能夠得到結(jié)果.可見3種算法的運算時間相當(dāng).

5 結(jié)論

在對辛Euler法、第一種數(shù)值算法、第二種數(shù)值算法進行擺方程實驗之前,通過理論的分析比較知道3個數(shù)值算法在算法的理論上的收斂性、局部穩(wěn)定性、誤差估計3個方面都是相當(dāng)?shù)?從算法的表面迭代格式上看,由Euler法和改進Euler法組合而成的第一種數(shù)值算法的迭代格式要比辛Euler法的迭代格式更好一些.首先,在算法實際迭代出算法格式的過程中,辛Euler法要比另外兩種算法更好,辛Euler法是構(gòu)造迭代格式最簡單的算法,其次是第一種數(shù)值算法,最后是第二種數(shù)值算法.盡管3種數(shù)值算法在理論上相當(dāng),但是在這3種數(shù)值算法的擺方程實驗中,可知辛Euler法具有很好的整體穩(wěn)定性,另外兩種數(shù)值算法的整體穩(wěn)定性不好,只是第一種數(shù)值算法的穩(wěn)定性比第二種數(shù)值算法的穩(wěn)定性要好一些.另外,3種數(shù)值算法的時間代價都是相當(dāng)?shù)?綜上所述,可知辛Euler法比另外兩種數(shù)值算法更好,更加具有實際應(yīng)用的價值.

[1] 余德浩,湯華中.微分方程數(shù)值解法[M].北京:科學(xué)出版社,2004.

[2] 李榮華,馮果忱.微分方程數(shù)值解法[M].3版.北京:高等教育出版社,1995.

[3] 戴嘉尊,邱建賢.微分方程數(shù)值解法[M].南京:東南大學(xué)出版社,2003.

[4] 毛濤濤,王正林,王玲.精通MATLAB GU I設(shè)計[M].北京:電子工業(yè)出版社,2008.

A New Algorithm for Second Order D ifferential Equations

TAN Jia-bo,CA ICheng-wei
(College of Info rm ation,Beijing W uzi University,B eijing101149)

Finds another two similar algorithm formats while constructing symplectic Euler method format through studying the process of symplectic Euler numerical algorithm.Comparing these two for mats with symplectic Euler method,finds the standard of an excellent numerical algorithm and the approach of improving such excellent numerical algorithm.By pendulum equation compares the numerical algorithms.Proves that the symplectic Euler algorithm format is better than the other two similar algorithm for mats.

second order differential equation;time cost;symplectic Eulermethod;new numerical algorithm;pendulum equation

O175.1

A

1007-0834(2011)01-0016-08

10.3969/j.issn.1007-0834.2011.01.007

2010-04-13

北京市優(yōu)秀教學(xué)團隊——“數(shù)學(xué)公共基礎(chǔ)系列課程教學(xué)團隊”支持項目(PHR200907230)

譚加博(1982—),男,黑龍江雞西人,北京物資學(xué)院信息學(xué)院講師,理學(xué)博士.