劉石威，劉文

(1.棗莊學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，山東 棗莊 277160；2.棗莊市 第四十一中學(xué)，山東 棗莊 277100)

0 引言

下面格式的二階微分方程問題一直廣受人們關(guān)注：

(1)

尤其當(dāng)它的解具有一定的周期性或振蕩性時.在化學(xué)物理學(xué)、物理化學(xué)、電子學(xué)、量子力學(xué)、天體力學(xué)等不同的領(lǐng)域中這類問題經(jīng)常出現(xiàn)，如果不能求出這類問題的解析解，那么求它的數(shù)值解顯得尤為重要.而三角擬合是其中較常使用的方法.Gautschi[1]與Lyche[2]首先給出此項技術(shù)的理論基礎(chǔ)，后來Raptis和Allison[3]針對二階微分方程構(gòu)造了一個Numerov類型的指數(shù)擬合方法，新方法對于解Schr?dinger類型的方程具有更高的效率.

1 三角擬合St?rmer-Verlet二步方法

對于求數(shù)值解的問題，我們考慮下面形式的顯示兩步格式：

yn+1=α1yn+α2yn-1+h2(α3f(xn,yn)+α4f(xn-1,yn-1)),

(2)

這個格式可以在[4]中找到.

讓方法(2)精確積分下面函數(shù)的線性組合:

{cos(ωx),sin(ωx)}

(3)

可得如下方程組：

(4)

在這里u=ωh.

α4=0.

(5)

上面系數(shù)的泰勒展式為:

α4=0.

(6)

以(5)所給系數(shù)的三角擬合St?rmer-Verlet兩步方法記為ST2F1，該方法的局部截斷誤差為:

(7)

因此這個方法的代數(shù)階為2.根據(jù)(6)，因為u=ωh,可知當(dāng)ω→0時，所得到的新方法ST2F1將變?yōu)樵嫉腟t?rmer-Verlet兩步方法.

2 相屬性與穩(wěn)定性分析

將得到的新方法應(yīng)用到試驗方程：

y″=λy,

(8)

可得到差分方程：

yn+1-2S(H2,u)yn+yn-1=0,

(9)

其中，S(H2,u)=(2-α3H2)/2,

H=λh，α3由(5)給出.

定義1:對于含S(H2,u)的新方法ST2F1，稱

φ(H2,u)=H-arccos(S(H2,u)),

為相滯.

稱ST2F1的相滯階數(shù)為q，當(dāng)

φ(H2,u)=O(Hq+1).

作替換k=H/u,k∈□,k≠1,此時這個新方法的相滯可以表示為

因此，這個新方法的相滯階數(shù)為2.

定義2:稱一個在平面H-u內(nèi)的二維平面區(qū)域D為含有S(H2,u)的方法ST2F1的絕對穩(wěn)定性區(qū)域，當(dāng)對于所有的(H,u)∈D，滿足

|S(H2,u)|<1,

并稱滿足

|S(H2,u)|=1.

的封閉曲線為絕對穩(wěn)定區(qū)域邊界.

在圖1中，我們給出新方法ST2F1的絕對穩(wěn)定性區(qū)域.

圖1 新方法ST2F1的絕對穩(wěn)定性區(qū)域

3 數(shù)值實驗

本節(jié)將新方法應(yīng)用于解四個振蕩問題，前三個是非齊次方程，最后一個是雙頻類周期問題.

選用下列四個方法予以比較:

原始的St?rmer-Verlet二步方法，見[4]，用ST2表示.

三角擬合St?rmer-Verlet二步方法，即本文所得到的新方法，用ST2F1表示.

原始的St?rmer類型的三步方法，見[4]，用ST3表示.

FSAL屬性的RKN方法，見[5]，用FRKN表示.

3.1 非齊次方程

求解下面問題:

y″=-100y+99sin(x),
y(0)=1,y′(0)=11.

(10)

其精確解為y(x)=cos(10x)+sin(10x)+sin(x).

對方程(10)做數(shù)值實驗，我們選擇區(qū)間0≤x≤100,針對于新方法ST2F1和FRKN選擇ω=10.評價的標準是比較四個方法的全局誤差及計算所用的步長，在圖2中我們給出數(shù)值演示結(jié)果.

3.2 非齊次方程

求解下面問題:

(11)

對方程(11)做數(shù)值實驗，我們選擇區(qū)間0≤x≤100,針對于新方法ST2F1和FRKN選擇ω=9.評價的標準是比較四個方法的全局誤差及計算所用的步長，在圖3中我們給出數(shù)值演示結(jié)果.

3.3 非齊次方程

求解下面問題:

(12)

對方程(12)做數(shù)值實驗，我們選擇區(qū)間0≤x≤100,針對于新方法ST2F1和FRKN選擇ω=13.評價的標準是比較四個方法的全局誤差及計算所用的步長，在圖4中我們給出數(shù)值演示結(jié)果.

3.4 雙頻類周期問題

求解下面問題:

y″=-169y+(480-160i)e3ix,
y(0)=4+i,y′(0)=-23+22i.

(13)

它等價于

u''(x)+169u(x)-480cos(3x)-160sin(3x)=0,

v''(x)+169v(x)-480sin(3x)+160cos(3x)=0.

其中y(x)=u(x)+iv(x).

其精確解為

u(x)=-2sin(13x)+cos(13x)+3cos(3x)+sin(3x),

v(x)=sin(13x)+2cos(13x)+3sin(3x)-cos(3x).

對方程(13)做數(shù)值實驗，我們選擇區(qū)間0≤x≤100,針對于新方法ST2F1和FRKN選擇ω=13.評價的標準是比較四個方法的全局誤差及計算所用的步長，在圖5中我們給出數(shù)值演示結(jié)果.

圖2 問題4.1數(shù)值演示結(jié)果 圖3 問題4.2數(shù)值演示結(jié)果

Fig.2EfficiencycurvesinProblem4.1Fig.3EfficiencycurvesinProblem4.2

圖4 問題4.3數(shù)值演示結(jié)果 圖5 問題4.4數(shù)值演示結(jié)果

Fig.4EfficiencycurvesinProblem4.3Fig.5EfficiencycurvesinProblem4.4

5 結(jié)論

從數(shù)值實驗的結(jié)果可知，三角擬合St?rmer-Verlet顯式二步方法在求解帶初值的振蕩問題時較其它三個方法更為高效.這是因為這個三角擬合方法能精確積分{cos(ωx),sin(ωx)}的線性組合，而其它的方法不具有這個特性.

參考文獻

[1]Gautschi W. ,Numerical integration of ordinary differential equations based on trigonometric polynomials[J],NumerischeMathematik, 1961:381-397.

[2]Lyche T.,Chebyshevian multistep methods for ordinary differential equations[J].NumerischeMathematik, 1972:65-75.

[3]Raptis A., Allison A.,Exponential-fitting methods for the numerical solution of the Schrodinger equation[J].ComputerPhysicsCommunications,1978: 1-5.

[4]Lambert J.D, Computational methods in ordinary differential equations[M].John Wiley & Sons Inc, New York, 1973.

[5]Van de Vyver H.A Runge-Kutta-Nystr?m pair for the numerical integration of perturbed oscillators, Computer Physics Communications, 2005:129-142.