秦于越 鄧子辰 胡偉鵬

(1．西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院，西安 710072)(2．大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116023)

諧振子的辛歐拉分析方法*

秦于越1?鄧子辰1，2胡偉鵬1

(1．西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院，西安 710072)(2．大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116023)

針對理想簡諧振子力學(xué)模型，研究了其守恒律，并利用辛歐拉格式分析簡諧振子振動過程．首先給出了諧振子系統(tǒng)的平方守恒律、周期守恒律和相差守恒律．構(gòu)造了諧振子的普通歐拉格式和辛歐拉格式，研究了兩種格式下三種守恒律各自的保持情況．模擬結(jié)果顯示:辛歐拉格式能夠精確保持時域守恒律(平方守恒律)，但無法保持頻域守恒律(周期守恒律和相差守恒律)．如要克服辛歐拉格式的不足，需按邢譽峰教授提出的方法進(jìn)行校正．

哈密爾頓， 保結(jié)構(gòu)算法， 辛歐拉， 簡諧振子， 守恒律

引言

簡諧振子是固體力學(xué)領(lǐng)域最簡單的理想力學(xué)模型之一．由于忽略系統(tǒng)的阻尼等耗散因素，簡諧振子的振動過程是一個典型的保守系統(tǒng)，振動能量應(yīng)該是守恒的，這一點在理論上早已得到了學(xué)術(shù)界的公認(rèn)，但是從數(shù)值分析中研究這類問題具有更廣闊的工程應(yīng)用背景．

1984年，保結(jié)構(gòu)幾何積分算法創(chuàng)始人在雙微國際會議上作了題為“On difference schemes and symplectic geometry”的大會報告［1］，首次為物理力學(xué)領(lǐng)域的數(shù)值算法提出了在數(shù)值求解過程中應(yīng)該充分關(guān)注物理力學(xué)系統(tǒng)的固有幾何性質(zhì)的保持，得到了國際同行的廣泛認(rèn)同，馮康先生也因此獲得了1997年自然科學(xué)一等獎．此后，辛算法在諸多領(lǐng)域得到了應(yīng)用:鐘萬勰院士將辛算法引入應(yīng)用力學(xué)領(lǐng)域的研究，并創(chuàng)造性地提出了精細(xì)積分算法用于時間歷程的積分運算［2］;秦孟兆等人繼承了馮康先生辛算法思想，提出了多步辛算法的構(gòu)造方法，并將其應(yīng)用于波動問題的數(shù)值分析過程，得到了一些有意義的結(jié)果［3，4］;Hairer教授等人全面總結(jié)辛算法研究成果，系統(tǒng)提出了針對常微分方程的保結(jié)構(gòu)幾何積分算法理論體系［5］;Budd 和 Piggott在 Hamilton框架下，系統(tǒng)闡述了保結(jié)構(gòu)算法的基本思想，為后續(xù)保結(jié)構(gòu)算法的發(fā)展指明了方向［6］;正是受保結(jié)構(gòu)算法思想的啟發(fā)，Bridges教授等人將針對有限維Hamilton動力學(xué)系統(tǒng)的辛算法推廣至無限維，創(chuàng)立了多辛算法［7］;鄧子辰教授等人針對非保守系統(tǒng)給出了廣義多辛算法［8－11］．

本文針對簡諧振子模型，研究了簡諧振子的固有幾何性質(zhì)及其對應(yīng)的守恒律，構(gòu)造辛歐拉格式模擬簡諧振子的振動過程，通過與常規(guī)歐拉格式的模擬結(jié)果進(jìn)行對比，充分展現(xiàn)辛歐拉格式的保結(jié)構(gòu)特點和長時間數(shù)值穩(wěn)定性特點．

1 簡諧振子的守恒律

考慮以下簡諧振子諧振問題

這是一個極其簡單的數(shù)學(xué)物理問題，其精確解可以很容易地得到，其中一個最為簡單的精確解為:

這個精確解很顯然是周期有界的，結(jié)合代數(shù)理論，得到其滿足的守恒律如下:

①平方守恒律:

這一守恒律表明諧振子的振動應(yīng)該始終都在一個閉合的圓上．

②周期守恒律:

其中T(u)和T(v)分別表示兩個廣義位移分量的振動周期．這一守恒律表明諧振子的兩個廣義位移分量的振動周期均是2π;

③相差守恒律:

需要說明的是以上守恒律是不依賴時間的，也就是說，按照保結(jié)構(gòu)思想，經(jīng)過無限長的時間后，以上三個守恒律仍然是應(yīng)該得到精確滿足的．

2 諧振問題的辛歐拉離散格式

令時間步長為h，依照差分理論，諧振問題(1)的普通歐拉差分離散格式為:

與嚴(yán)格保能量的歐拉中點辛差分［12］格式相比，格式(6)是一個顯式格式，計算速度會更快．從這一差分格式中很容易就可以得到第n+1步的結(jié)果與第n步的結(jié)果之間的代數(shù)關(guān)系為:

這一代數(shù)關(guān)系不僅不依賴于時間，而且不依賴于初始條件和邊界條件．代數(shù)關(guān)系(7)從根本上與平方守恒律(3)相背．同時，無論時間步長多么小，由普通歐拉格式(6)得到的諧振子的振動將是無界的，即，隨著時間的演化，振動幅值趨近于無窮大，這從根本上破壞了原系統(tǒng)的周期性，更談不上數(shù)值解能滿足周期守恒律(4)和相差守恒律(5)了．

注意到系統(tǒng)(1)是一個典型的分離形式的Hamilton系統(tǒng)，其Hamilton數(shù)為:

采用辛歐拉離散方法［5］，離散系統(tǒng)(1)，得到離散辛格式:

其迭代雅克比矩陣為:

容易驗證該雅克比矩陣是辛的．

3 數(shù)值試驗

假定諧振子的初始條件為:

圖1 諧振子兩個廣義位移分量變化情況Fig．1 The evolutions of functions and

取時間步長h=0．05s，分別采用普通歐拉格式(6)和辛歐拉格式(9)模擬諧振子的振動過程，得到兩個廣義位移分量結(jié)果如圖1，同時得到三個守恒律的保持情況如圖2和圖3(后兩種守恒律的保持情況只能給出辛歐拉格式的結(jié)果，因為普通歐拉格式得到的結(jié)果是非周期的)．

圖2給出了兩種格式得到的數(shù)值結(jié)果在平方守恒律方面的結(jié)果對比，從對比結(jié)果不難發(fā)現(xiàn):普通歐拉格式得到的結(jié)果滿足方程(7)，即u2+v2的數(shù)值線性增長;而辛歐拉格式的u2+v2數(shù)值始終保持為1．

圖3 周期守恒律和相差守恒律結(jié)果Fig．3 The conservation law of period and phase difference

從周期守恒律和相位差守恒律結(jié)果可以看出:隨著模擬的進(jìn)行，諧振子兩個廣義位移分量的周期緩慢線性增長，同時兩個廣義位移之間的相位差也緩慢線性增長，這說明辛算法不能精確保持系統(tǒng)的頻域特性．這一點與邢譽峰教授關(guān)于辛算法的相位誤差的研究結(jié)果［13］吻合，如果要精確保持系統(tǒng)的頻域特性，需要按照邢譽峰教授建議的方法進(jìn)行校正．

4 結(jié)論

本文針對諧振子理想模型，分析其多種守恒律，并從理論和數(shù)值模擬結(jié)果兩個方面分別比較了諧振子模型普通歐拉格式和辛歐拉格式在保持守恒律方面的差異．對比結(jié)果表明:辛歐拉格式能夠較好的長時間模擬諧振子的振動過程，同時長時間精確保持諧振子時域內(nèi)的平方守恒律，但是不能精確保持頻域內(nèi)的周期守恒律和相位差守恒律，若要求其能夠精確保持頻域內(nèi)的相關(guān)守恒律，還需要對辛歐拉格式進(jìn)行校正．

1 Feng K．On difference schemes and symplectic geometry，Proceeding of the 1984 Beijing symposium on D．D．，Beijing:Science Press，1984，42～58

2 Zhong W X．Some developments of computational solid Mechanics in China．Computers＆Structures，1988，30(4):783～788

3 Qin M Z，Zhang M Q．Multi－stage symplectic schemes of two kinds of Hamiltonian systems for wave equations．Computers＆Mathematics with Applications，1990，19(10):51～62

4 Zhao P F，Qin M Z．Multisymplectic geometry and multisymplectic Preissmann scheme for the KdV equation．Journal of Physics a－Mathematical and General，2000，33(18):3613～3626

5 Hairer E，Lubich C，Wanner G．Geometric numerical integration:structure preserving algorithms for Ordinary Differential Equations．Berlin:Springer－Verlag，2002

6 Budd C J，Piggott M D，Geometric integration and its applications．Amsterdam:Handbook of Numerical Analysis，2003:35～139

7 Bridges T J．Multi－symplectic structures and wave propagation．Mathematical proceedings of the Cambridge Philosophical Society，1997，121(1):147～190

8 Hu W P，Deng Z C，Han S M，Zhang W R．Generalized multi－symplectic integrators for a class of Hamiltonian nonlinear wave PDEs．Journal of Computational Physics，2013，235:394～406

9 胡偉鵬，鄧子辰．大阻尼桿振動的廣義多辛算法．動力學(xué)與控制學(xué)報，2013，11(1):1～4(Hu W P，Deng Z C．Generalized multi－symplectic method for vibration of big damping bar．Journal of Dynamics and Control，2013，11(1):1～4(in Chinese))

10 Hu W P，Deng Z C，Wang B，Ouyang H J．Chaos in an embedded single－walled carbon nanotube．Nonlinear Dynamics，2013，72(1 －2):389 ～398

11 Hu W P，Deng Z C，Ouyang H J．Generalized Multisymplectic Method for Dynamic Responses of Continuous Beam under Moving Load．International Journal of Applied Mechanics，2013，5(3):1350033

12 邢譽峰，楊蓉．動力學(xué)平衡方程的Euler中點辛差分求解格式．力學(xué)學(xué)報，2007，39(1):100～105(Xing Y F，Yang R．Application of Euler midpoint symplectic integration method for the solution of dynamic equilibrium equations．Acta Mechanica Sinica，2007，39(5):100～105(in Chinese))

13 邢譽峰，楊蓉．單步辛算法的相位誤差分析及修正．力學(xué)學(xué)報，2007，39(5):668～671(Xing Y F，Yang R．Phase errors and their correction in symplectic implicit sin－gle－step algorithm．Acta Mechanica Sinica，2007，39(5): 668 ～671(in Chinese))

*The project supported by the National Natural Science Foundation of China(11172239，11002115，11372253)，the Doctoral Program Foundation of Education Ministry of China(20126102110023)and the Open Foundation of State Key Laboratory of Structural Analysis of Industrial Equipment(GZ0802)

? Corresponding author E－mail:769482448@qq．com

SYMPLECTIC EULER METHOD FOR HARMONIC OSCILLATOR*

Qin Yuyue1?Deng Zichen1，2Hu Weipeng1
(1Department of Engineering Mechanics，Northwestern Polytechnical University，Xi’an710072，China)(2State Key Laboratory of Structural Analysis of Industrial Equipment，Dalian University of Technology，Dalian116023，China)

Focusing on the conservation properties，the symplectic Euler scheme of the harmonic oscillator was constructed to analyze its vibration properties．Firstly，three conservation laws，including the square conservation law，the period conservation law and the phase difference conservation law，were presented for the harmonic oscillator．And then，the common Euler scheme and the symplectic Euler scheme were constructed to study the above three conservation laws．The numerical results imply that the symplectic Euler scheme can preserve the conservation law in time domain(the square conservation law)exactly，but can’t preserve the conservation laws in phase domain(the period conservation law and the phase difference conservation law)，which is the shortcoming of the symplectic method but can be overcome by the modification method presented by Prof．Xing．

Hamiltonians， structure － preserving method， symplectic Euler， harmonic oscillator， conservation law

21 January 2013，

29 August 2013．

10．6052/1672－6553－2013－094

2013－01－21 收到第 1 稿，2013－08－29 收到修改稿．

*國家自然科學(xué)基金(11172239，11002115，11372253)，博士點基金(20126102110023)和大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室開放基金(GZ0802)

E－mail:769482448@qq．com