求解無(wú)阻尼Duffing方程的三角擬合顯式對(duì)稱六步法

房永磊,郭娟

(棗莊學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山東棗莊　277100)

?

求解無(wú)阻尼Duffing方程的三角擬合顯式對(duì)稱六步法

房永磊,郭娟

(棗莊學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山東棗莊277100)

本文利用多頻三角擬合多步法求解無(wú)阻尼Duffing方程，構(gòu)造包含不同共振譜的四種三角擬合對(duì)稱六步法，分析新方法的穩(wěn)定性和相性質(zhì).數(shù)值實(shí)驗(yàn)的結(jié)果顯示新方法的有效性.

顯式對(duì)稱方法；多頻問(wèn)題；相分析①

0　引言

無(wú)阻尼Duffing方程在經(jīng)典力學(xué)中具有非常重要的應(yīng)用.無(wú)阻尼Duffing方程的一般形式為：

(1)

這個(gè)方程沒(méi)有精確解，Micken表明等式(1)的解可以擴(kuò)展為系列周期函數(shù)形如[1]：

(2)

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

然而, 估計(jì)出所有系數(shù)A2i+1的值是一個(gè)很大的挑戰(zhàn).Duffing方程的數(shù)值積分問(wèn)題的研究越來(lái)越多.最近，Wang[2]通過(guò)應(yīng)用傅里葉頻譜的高階導(dǎo)數(shù)，構(gòu)造一個(gè)新的三角擬合Obrechkoff單步法，這個(gè)方法是隱式的，但每一個(gè)積分步長(zhǎng)都需要求解一個(gè)非線性方程組.這對(duì)于無(wú)阻尼Duffing方程來(lái)說(shuō)，計(jì)算成本很高.

本文的研究目標(biāo)是針對(duì)無(wú)阻尼Duffing方程給出一類三角擬合顯式對(duì)稱六步法.結(jié)構(gòu)如下：第二節(jié)通過(guò)傅立葉頻的擬合技術(shù)，構(gòu)造了三個(gè)新的三角擬合對(duì)稱六步法；第三節(jié)介紹了新方法對(duì)于非阻尼Duffing方程的局部截?cái)嗾`差；第四節(jié)分析了新方法的穩(wěn)定性和相性質(zhì)；第五節(jié)給出了數(shù)值實(shí)驗(yàn)顯示新方法的有效性；第六節(jié)給出了一些結(jié)論.

1　新方法的構(gòu)造

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

(3)

考慮求解問(wèn)題(3)的顯式對(duì)稱六步法：

(4)

這個(gè)方法代數(shù)階是6，把它定義Method I.為推導(dǎo)出適應(yīng)振蕩問(wèn)題的新方法，可以讓方法(4)精確積分以下線性組合

(5)

1.1方法二的構(gòu)造

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

(6)

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

該方法的局部截?cái)嗾`差為：

這個(gè)方法代數(shù)階是6，記為Method II.

1.2方法三的構(gòu)造

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

(7)

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

這個(gè)方法的代數(shù)階是6，記為Method III.

1.3方法四的構(gòu)造

讓方法(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)

這個(gè)方法代數(shù)階是6，記為Method IV.

2　誤差分析

本節(jié)針對(duì)問(wèn)題(1)將給出方法I-IV的局部截?cái)嗾`差，取參數(shù)B=0.002，w=1.01，初值為：y(0)=0.200426728069666,y′(0)=0,問(wèn)題(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.

以函數(shù)(9)作為非阻尼Duffing方程的精確解，方法I-IV的局部截?cái)嗾`差可記為：

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

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

3　穩(wěn)定性分析

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

3.1對(duì)稱多步法的穩(wěn)定性

把對(duì)稱步方法

(10)

應(yīng)用到實(shí)驗(yàn)方程

y''=-λ2y

(11)

得到下列差分方程

(12)

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

(13)

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

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

定理1[10]： 稱

(14)

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

4.2新方法的穩(wěn)定性

把方法II-IV應(yīng)用到特征方程可以得到新方法的特征方程為：

(15)

為此有以下定義.

定義3[9]：稱平面上的區(qū)域?yàn)橐?15)為 特征方程的新方法的穩(wěn)定性區(qū)域，如果滿足

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

為新方法的相延遲誤差.

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

圖1　方法II的穩(wěn)定性區(qū)域

4　數(shù)值試驗(yàn)

從表1中我們可以看到，新方法II，III和IV對(duì)于無(wú)阻尼Duffing方程式非常有效的.方法III和IV比TFNT方法更加準(zhǔn)確，尤其當(dāng)步長(zhǎng)時(shí)，方法Ⅳ仍保持較高精度，而其它方法都出現(xiàn)了較大誤差.

圖2　方法III的穩(wěn)定性區(qū)域

圖3　方法IV的穩(wěn)定性區(qū)域

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　結(jié)論

針對(duì)無(wú)阻尼Duffing方程，通過(guò)采用傅里葉譜擬合前幾個(gè)分量,本文給出一類新的顯式三角擬合對(duì)稱六步法，分析了新方法的穩(wěn)定性和相位特性，數(shù)值實(shí)驗(yàn)的結(jié)果驗(yàn)證了新方法的高效.方法IV高效性說(shuō)明它已經(jīng)很好的擬合了Duffing方程傅里葉頻譜的前三個(gè)分量.如果考慮更多傅里葉譜分量擬合可以構(gòu)建更多精確的三角擬合積分方法.

[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.

[責(zé)任編輯：呂海玲]

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

國(guó)家自然科學(xué)基金(項(xiàng)目編號(hào)：11571302).

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

O241.82

A

1004-7077(2016)05-0009-06