求解無阻尼Duffing方程的三角擬合顯式對稱六步法

房永磊,郭娟

(棗莊學院數學與統計學院，山東棗莊　277100)

?

求解無阻尼Duffing方程的三角擬合顯式對稱六步法

房永磊,郭娟

(棗莊學院數學與統計學院，山東棗莊277100)

本文利用多頻三角擬合多步法求解無阻尼Duffing方程，構造包含不同共振譜的四種三角擬合對稱六步法，分析新方法的穩定性和相性質.數值實驗的結果顯示新方法的有效性.

顯式對稱方法；多頻問題；相分析①

0　引言

無阻尼Duffing方程在經典力學中具有非常重要的應用.無阻尼Duffing方程的一般形式為：

(1)

這個方程沒有精確解，Micken表明等式(1)的解可以擴展為系列周期函數形如[1]：

(2)

從等式(1)的通解展開式(2)看出，等式(1)的頻譜都是奇次諧頻ω,3ω,···,共振現象便是由這些諧頻產生.

然而, 估計出所有系數A2i+1的值是一個很大的挑戰.Duffing方程的數值積分問題的研究越來越多.最近，Wang[2]通過應用傅里葉頻譜的高階導數，構造一個新的三角擬合Obrechkoff單步法，這個方法是隱式的，但每一個積分步長都需要求解一個非線性方程組.這對于無阻尼Duffing方程來說，計算成本很高.

本文的研究目標是針對無阻尼Duffing方程給出一類三角擬合顯式對稱六步法.結構如下：第二節通過傅立葉頻的擬合技術，構造了三個新的三角擬合對稱六步法；第三節介紹了新方法對于非阻尼Duffing方程的局部截斷誤差；第四節分析了新方法的穩定性和相性質；第五節給出了數值實驗顯示新方法的有效性；第六節給出了一些結論.

1　新方法的構造

Duffing方程(1)可以寫為二階常微分方程的一般形式[3]：

(3)

考慮求解問題(3)的顯式對稱六步法：

(4)

這個方法代數階是6，把它定義Method I.為推導出適應振蕩問題的新方法，可以讓方法(4)精確積分以下線性組合

(5)

1.1方法二的構造

按照文獻[4,5,6,7]，要求方法(4)精確積分的值，可得到：

(6)

取b1=-1/6,b2=61/24, 解(6)得到

該方法的局部截斷誤差為：

這個方法代數階是6，記為Method II.

1.2方法三的構造

要求方法(4)精確積分，可得到：

(7)

取b2=61/24, 求解方程組(7)得到：

這個方法的代數階是6，記為Method III.

1.3方法四的構造

讓方法(4)精確積分，可得到：

(8)

求解方程組(8)，可得到：

b0=(csc(u)6sec(u)3(-225cos(3u)2+25cos(5u)(cos(6u)-cos(9u))

+9cos(3u)(25cos(25u)-cos(10u)+cos(15u))+cos(u)(-25cos(6u)

+25cos(9u)+9cos(10u)-9(15u)+1800cos(2u)2sin(u)2)))/N

b1=((225cos(3u)(cos(6u)-cos(10u))-16cos(6u)cos(10u)+9cos(10u)

+cos(2u)(-200cos(6u)-25cos(9u)+216cos(10u)+9cos(15u))

+25cos(9u)cos(10u)-9cos(6u)cos(15u)csc(u)6sec(u)3)/N

b2=((2432+4504cos(u)+3559cos(2u)+2434cos(3u)+1579cos(4u)+1064cos(5u)+829cos(6u)+734cos(7u)+584cos(8u)+354cos(9u)+164cos(10u)+54cos(11u)

這個方法代數階是6，記為Method IV.

2　誤差分析

本節針對問題(1)將給出方法I-IV的局部截斷誤差，取參數B=0.002，w=1.01，初值為：y(0)=0.200426728069666,y′(0)=0,問題(1)的精確的近似解:

(9)

其中， A1=0.2001794775366180,A3=0.000246946143255837,A5=3.04014985249×10-7,A7=3.744×10-10,A9=3.74349084×10-10,A11=5.68×10-16.

以函數(9)作為非阻尼Duffing方程的精確解，方法I-IV的局部截斷誤差可記為：

由A1≈0.2,A3≈10-3,···,A11≈10-16,可以得到

由以上可以得出,通過傅里葉三角擬合技術得到的數值方法具有很好的誤差結果.

3　穩定性分析

本節將分析上文提到新方法的穩定性和相對特性.Lambert 和Watson的穩定理論[8]被應用到Colleman和Ixaru提出的求解y''=f(x,y)指數擬合對稱方法[9].

3.1對稱多步法的穩定性

把對稱步方法

(10)

應用到實驗方程

y''=-λ2y

(11)

得到下列差分方程

(12)

其中， θ=λh. 方法(10)的特征多項式可記為：

(13)

定義1[8]：稱以(12)為特征方程的步方法(10)的周期區間為，如果差分方程的特征根滿足

其中a(θ)是θ的實值函數.

定理1[10]： 稱

(14)

為以(12)為特征方程的2k步方法(10)的相誤差.

4.2新方法的穩定性

把方法II-IV應用到特征方程可以得到新方法的特征方程為：

(15)

為此有以下定義.

定義3[9]：稱平面上的區域為以(15)為 特征方程的新方法的穩定性區域，如果滿足

其中a(u,θ)是u,θ的實值函數.

為新方法的相延遲誤差.

因此，方法II-IV相位延遲都為6階.圖1-3給出方法II-IV的穩定性區域，從圖中可以看出新方法II-IV的穩定性區域依次增大.

圖1　方法II的穩定性區域

4　數值試驗

從表1中我們可以看到，新方法II，III和IV對于無阻尼Duffing方程式非常有效的.方法III和IV比TFNT方法更加準確，尤其當步長時，方法Ⅳ仍保持較高精度，而其它方法都出現了較大誤差.

圖2　方法III的穩定性區域

圖3　方法IV的穩定性區域

hTFNTMethodIMethodIIMethodIIIMethodIV2NaNNaNNaNNaN2.9e-810.15NaNNaNNaN1.4e-70.58.7e-59.6e-41.7e-61.4e-51.1e-100.251.3e-68.6e-62.2e-68.5e-72.3e-100.1251.9e-81.3e-73.1e-81.3e-91.6e-10

5　結論

針對無阻尼Duffing方程，通過采用傅里葉譜擬合前幾個分量,本文給出一類新的顯式三角擬合對稱六步法，分析了新方法的穩定性和相位特性，數值實驗的結果驗證了新方法的高效.方法IV高效性說明它已經很好的擬合了Duffing方程傅里葉頻譜的前三個分量.如果考慮更多傅里葉譜分量擬合可以構建更多精確的三角擬合積分方法.

[1]R.E. Mickens, An Introduction to Nonlinear Oscillations[M]. New York, Cambridge University Press,1981.

[2]Z.C. Wang, Obrechkoff one-step method fitted with Fourier spectrum for undamped Duffing equation[J]. Computer Physics Communication,2006 ,175 (11-12): 692-699.

[3]T.E. Simos, Exponentially-fitted and trigonometrically-fitted methods for long-term integration of orbital problems[J]. New Astronomy. 2002,7 (1): 1-7.

[4]G. Vanden Berghe, H. De Meyer, M. Van Daele, T. Van Hecke,Exponentially fitted Runge-Kutta methods[J].Computer Physics Communication, 1999, 123 (1-2): 7-15.

[5]J.M. Franco, Exponentially fitted explicit Runge-Kutta-Nystrm methods[J]. Journal of Computational and Applied Mathematics,2004, 167 (1):1-19.

[6]Y.L. Fang, X.Y. Wu, A trigonometrically fitted explicit Numerov-type method for second-order initial value problems with oscillating solutions[J].Applied Numerical Mathematics, 2008,58 (3) 341-351.

[7]Z.C. Wang, Trigonometrically-fitted method for a periodic initial valueproblem with two frequencies[J].Computer Physics Communication, 2006,175 (4) :241-249.

[8]J.D. Lambert, I.A. Watson, Symmetric multistep methods for periodic initial value problems[J].Journal of the Institute of Mathematics and its Applications,1976,18 (2): 189-202.

[9]J.P. Coleman, L.Gr. Ixaru, P-stability and exponential-fitting methods fory''=f(x,y) [J].IMA Journal of Numerical Anaysis,1996 16 (2) :179-199.

[10]G.A. Panopoulos, Z.A. Anastassi, T.E. Simos, A symmetric eight-step predictor-corrector method for the numerical solution of the radial Schringer equation and related IVPs with oscillating solutions[J]. Computer Physics Communication. 2011,82 (8): 1626-1637.

[11]Y.L. Fang, Y.Z. Song, X.Y. Wu, Trigonometrically fitted explicit Numerov-type method for periodic IVPs with two frequencies[J]. Computer Physics Communication, 2008,179 (11) : 801-811.

[責任編輯：呂海玲]

Trigonometrically Fitted Explicit Symmetric Six-Step Methods for Undamped Duffing Equation

FANG Yong-lei,GUO Juan

(School of Mathematics and Statistics, Zaozhuang University, Zaozhuang 277160,China)

In this paper, a class of trigonometrically fitted explicit symmetric six-step methods for the numerical integration of undamped duffing equation is developed. Using multi-frequency trigonometrically fitted technique, we design four types of six-step methods which contain different resonance spectrum. Numerical stability and phase properties of the new methods are analyzed to explain the final numerical experiment.

explicit symmetric method; multi-frequency problems; undamped duffing equation

2016-08-19

國家自然科學基金(項目編號：11571302).

房永磊(1981-)，山東棗莊人，棗莊學院數學與統計學院教授，理學博士，主要從事微分方程數值解的研究.

O241.82

A

1004-7077(2016)05-0009-06