求解无阻尼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