非線性常微分方程高階諧波平衡法傅里葉展開的簡(jiǎn)化

唐元璋,翁雪濤，樓京俊，林雄偉，張暉

（1.海軍工程大學(xué)動(dòng)力工程學(xué)院，武漢430033；2.船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室，武漢430033）

非線性常微分方程高階諧波平衡法傅里葉展開的簡(jiǎn)化

唐元璋1,2,翁雪濤1,2，樓京俊1,2，林雄偉1,2，張暉1,2

（1.海軍工程大學(xué)動(dòng)力工程學(xué)院，武漢430033；2.船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室，武漢430033）

簡(jiǎn)化了一種求取非線性常微分方程高階諧波解的近似解析計(jì)算方法。對(duì)平方和立方非線性項(xiàng)的傅里葉展開過程進(jìn)行改進(jìn)和簡(jiǎn)化，使計(jì)算過程變?yōu)閮纱尉仃囘\(yùn)算即可完成展開過程，且兩次矩陣運(yùn)算過程一致，易于編程。以Duffing方程為算例，計(jì)算結(jié)果與數(shù)值方法一致，運(yùn)算效率有所提高。

振動(dòng)與波；非線性常微分方程；Duffing方程；傅立葉展開；諧波平衡法

線性隔振系統(tǒng)存在共振難避免等問題[1，2]。因此學(xué)者們對(duì)非線性隔振器開展了研究[3，4]，但與線性隔振系統(tǒng)相比，非線性隔振系統(tǒng)的響應(yīng)復(fù)雜得多。很多非線性振動(dòng)系統(tǒng)可以簡(jiǎn)化為單自由度質(zhì)點(diǎn)與非線性彈簧組合的振動(dòng)系統(tǒng)，這類系統(tǒng)的動(dòng)力學(xué)方程即為非線性常微分方程。圖1所示的系統(tǒng)動(dòng)力學(xué)方程無量綱化后稱為Duffing方程[5]。研究非線性常微分方程的主要方法包括解析方法和數(shù)值方法，數(shù)值方法只能得到解的時(shí)間序列，不能直接得到周期解各諧波分量的幅值與系統(tǒng)各參數(shù)之間的變化關(guān)系，所以解析方法仍然是重要的研究方法。非線性常微分方程的精確解析解大多無法求解，這里的解析方法主要是指近似解析方法，求解非線性常微分方程的近似解析方法有諧波平衡法、攝動(dòng)法、平均法、多尺度法和漸進(jìn)法等[6，7]。各種近似解析方法中，諧波平衡法概念最簡(jiǎn)單明了，且適用于強(qiáng)非線性系統(tǒng)，基本思想是將振動(dòng)系統(tǒng)的激勵(lì)項(xiàng)和方程的解都展開成傅立葉級(jí)數(shù)，為了保證系統(tǒng)的作用力和慣性力的各階諧波分量自相平衡，必須令動(dòng)力學(xué)方程兩端的同階諧波的系數(shù)相等，從而得到包含未知系數(shù)的一系列代數(shù)方程，以確定待定的傅立葉級(jí)數(shù)的系數(shù)。諧波平衡法解非線性常微分方程的過程就是確定待定的傅立葉級(jí)數(shù)的系數(shù)的過程。

圖1 非線性隔振系統(tǒng)

應(yīng)用諧波平衡法是假設(shè)非線性常微分方程的解只包含一個(gè)諧波[8，9]，本文簡(jiǎn)稱單一諧波解，這在研究共振時(shí)是有效的，但非共振情況下其他諧波項(xiàng)不能忽略，這需要研究高階諧波解，而對(duì)于高階諧波解的研究很少涉及[10，12]，而單一諧波或諧波個(gè)數(shù)較少，則無法求得高精度的解；如果非線性常微分方程包含平方項(xiàng)和立方項(xiàng)，計(jì)算高階諧波解的傅立葉級(jí)數(shù)系數(shù)時(shí)需將其展開，當(dāng)諧波階數(shù)超過5時(shí)，這個(gè)展開過程如果用直接積分的方法，將產(chǎn)生很大的計(jì)算量，計(jì)算機(jī)難以計(jì)算。為了減少積分運(yùn)算，Kawakami[13]等用求取偏導(dǎo)數(shù)的方法計(jì)算平方項(xiàng)和立方項(xiàng)傅里葉展開過程，但計(jì)算過程較為復(fù)雜；Ameer Hassan[14]用到一種直接展開計(jì)算的高階諧波平衡法，沈建和[15]等研究了類似的計(jì)算方法。這些方法，計(jì)算過程較為復(fù)雜，不利于編程實(shí)現(xiàn)。本文改進(jìn)了一種針對(duì)諧波平衡法傅里葉系數(shù)展開的計(jì)算方法，力求簡(jiǎn)化過程，使得計(jì)算平方項(xiàng)和立方項(xiàng)的展開過程統(tǒng)一，便于編程實(shí)現(xiàn)，也易于擴(kuò)展到計(jì)算其它類型非線性項(xiàng)的傅里葉展開。這一改進(jìn)的傅里葉展開計(jì)算過程不僅可以在諧波平衡法中應(yīng)用，也可以應(yīng)用到橢圓函數(shù)諧波平衡法[16]、增量諧波平衡法中[17]。為了檢驗(yàn)本文計(jì)算方法的正確性和有效性，對(duì)Duffing方程進(jìn)行了求解計(jì)算，并與數(shù)值方法做了對(duì)比。

1 計(jì)算過程

1.1 高階諧波平衡法的求解過程

第一步：將解設(shè)為傅立葉級(jí)數(shù)的解。設(shè)它的周期解為傅立葉級(jí)數(shù)形式

式(2)中Fc0，F(xiàn)ck，F(xiàn)sk是包含a0，ak，bk的多項(xiàng)式，其中k=1,2,3…∞。根據(jù)傅立葉級(jí)數(shù)系數(shù)的定義，定義P為求傅立葉級(jí)數(shù)系數(shù)的運(yùn)算，P與Fc0，F(xiàn)ck，F(xiàn)sk有如下關(guān)系

第三步：列出包含未知系數(shù)的代數(shù)方程組。由(2)式可得Fc0=0，F(xiàn)ck=0，F(xiàn)sk=0即

其中k=1,2,3…n,k取有限項(xiàng)即為近似解。這樣根據(jù)三角級(jí)數(shù)的正交性和式（1）—（8）可列出2n+1個(gè)方程。

第四步：解方程組。當(dāng)n很?。╪≤2）可求取幅值與頻率之間的解析關(guān)系，當(dāng)n很大時(shí)則只能用數(shù)值方法得到半數(shù)值半解析解，方程組（6）—（8）是2n+1個(gè)方程組成的，可用牛頓法或其他數(shù)值算法[18，19]解得各系數(shù)a0,a1…ak,b1…bk。

1.2 多項(xiàng)式函數(shù)的Fourier展開

為了列出1.1節(jié)第三步的代數(shù)方程組，如果直接按式(3)，(4)，(5)直接積分求解，那么計(jì)算量將非常大，n≥5時(shí)用Mathematica計(jì)算消耗內(nèi)存很大而無法運(yùn)算。如設(shè)更高階諧波則很難計(jì)算。為了減少計(jì)算量，本文對(duì)多項(xiàng)式函數(shù)Fourier展開法進(jìn)行了改進(jìn)，使得過程簡(jiǎn)單，易于理解和編程實(shí)現(xiàn)。設(shè)非線性常微分方程的近似解為有限項(xiàng)的傅立葉級(jí)數(shù)

n為式（9）中包含的諧波階數(shù)。c0,ck,dk是包含a0,ak,bk的多項(xiàng)式。利用式(6)，(7)，(8)可以求取c0,ck,dk與a0,ak,bk的關(guān)系。為了簡(jiǎn)化計(jì)算過程，改用矩陣計(jì)算，這樣也便于編程計(jì)算。設(shè)

c中共有4n+1項(xiàng)。則x(t)可表示為式（13）

定義矩陣

則有

根據(jù)式（17）可以很快得到x2(t)的傅立葉系數(shù)Pc0(x2(t))，Pck(x2(t))，Psk(x2(t))，也就是c0，ck，dk其中(k=1，2，3…2n)。

要求K(x2(t))，首先得研究TTT矩陣，計(jì)算這個(gè)矩陣等價(jià)于x2(t)的展開過程。

根據(jù)TTT的特點(diǎn)，利用三角級(jí)數(shù)正交性和積化和差公式即可計(jì)算K(x2(t))，用K(x2(t))[i，j]表示第i行j列。其中i,j=1，3…2n，各矩陣元素計(jì)算如下

當(dāng)i≠j時(shí)

其中Sign[]和Abs[]分別為取符號(hào)運(yùn)算和取絕對(duì)值運(yùn)算

類似地x3(t)的傅立葉系數(shù)P(x3(t))也可求出，但因?yàn)橹恍枇谐?n+1個(gè)方程，不必計(jì)算出全部P(x3(t))，所以另設(shè)

其中T1和c1只包含2n+1項(xiàng)

定義矩陣

則有

根據(jù)式（37）可以得到x3(t)的傅立葉系數(shù)Pc0(x3(t)),Pck(x3(t)),Psk(x3(t)),也就是e0,ek,fk,其中k=1,2,3…n。若要計(jì)算x3(t)全部的系數(shù)則將T1和c1設(shè)為6n+1項(xiàng)即可。

類似地可以得到以下結(jié)論

當(dāng)i≠j時(shí)

當(dāng)i=j時(shí)

綜上所述，Pc0(x3(t)),Pck(x3(t)),Psk(x3(t))即x3(t)的傅立葉系數(shù)e0,ek,fk可由K(x3(t)) c1得到e0,ek,fk與c0,ck,dk的關(guān)系；由K(x2(t)) c可得到c0,ck,dk與a0,ak,bk的關(guān)系。只需要通過兩次矩陣運(yùn)算就可就求得Pc0(x3(t)),Pck(x3(t)),Psk(x3(t))。這樣通過式(6)，(7)，(8)可以很容易得到2n+1個(gè)方程。

1.3 總結(jié)計(jì)算過程

整個(gè)計(jì)算過程總結(jié)為兩次矩陣運(yùn)算和解方程組運(yùn)算。

1.4 誤差

這里的誤差只討論近似解析解的計(jì)算誤差，因精確解尚無法求得，近似解析解與精確解的誤差對(duì)比不做討論。計(jì)算誤差來源于兩方面，一是牛頓法的計(jì)算誤差E1，可以通過增加迭代次數(shù)來減小誤差。

二是截?cái)嗾`差，所設(shè)諧波解僅為有限項(xiàng)n,多項(xiàng)式最高次項(xiàng)xm()t展開后最多有mn項(xiàng)，利用e0,ek, fk可求得截?cái)嗾`差E2。為了求出全部的系數(shù)只需將1.2節(jié)中T1和c1設(shè)為包含2mn+1項(xiàng)即可

總誤差為

2 算例

2.1 計(jì)算周期解

硬彈簧Duffing振子[1，20]是非線性振動(dòng)模型的典型代表，它的動(dòng)力學(xué)方程經(jīng)無量綱化后為

式（52）是一個(gè)多項(xiàng)式非線性方程，稱為Duffing方程，式中δ為阻尼率，α1和α3分別為線性項(xiàng)系數(shù)和立方非線性項(xiàng)系數(shù)，A為激勵(lì)力，ω激勵(lì)頻率,φ為相位；分別為加速度、速度、位移，以上各參數(shù)均為量綱為1的量。

聯(lián)合式（1）—（8）及式（52）很快可以得到方程組

固定一組參數(shù)進(jìn)行運(yùn)算，令α1=1，α3=1，A=1，φ=0，ω=0.9，δ=0.1。將諧波階數(shù)n設(shè)定為30進(jìn)行運(yùn)算,得到結(jié)果是高階諧波近似解析解的傅立葉系數(shù)，列于表1。

為了驗(yàn)證計(jì)算方法的正確性，與數(shù)值計(jì)算方法進(jìn)行對(duì)比，如圖2所示，實(shí)線表示數(shù)值計(jì)算的結(jié)果，實(shí)點(diǎn)表示本文方法，本文稱之為矩陣方法。兩種計(jì)算方法所得的結(jié)果是一致的。

2.2 誤差分析

按1.4節(jié)的定義計(jì)算誤差，各參數(shù)不變，只改變所設(shè)近似解析解的諧波階數(shù)。表2表明隨著諧波階數(shù)的增加，E1和E2都會(huì)下降，當(dāng)n=30時(shí)，E2＜＜E1這說明計(jì)算方法的精度是足夠的。

2.3 運(yùn)算效率分析

本文方法是將諧波平衡法的傅里葉展開過程整理成兩次矩陣運(yùn)算，非常適合在Matlab軟件中計(jì)算[21]，減少了代碼語句，如表3所示，運(yùn)算效率有一定的提高。

圖2 矩陣方法與數(shù)值方法的對(duì)比

表1 高階諧波解的傅立葉系數(shù)

表2 誤差分析

表3 運(yùn)算時(shí)間比較

3 結(jié)語

總結(jié)了多種近似解析方法后，針對(duì)帶有平方項(xiàng)和立方項(xiàng)的非線性常微分方程，簡(jiǎn)化和改進(jìn)了一種高階諧波平衡法。這種計(jì)算方法僅包含兩次矩陣運(yùn)算和一次解方程組運(yùn)算，兩次矩陣運(yùn)算的過程一致，步驟簡(jiǎn)單，可以方便地進(jìn)行編程計(jì)算。本文的主要工作包括：

（1）與多種近似解析方法相比，為了求得非線性常微分方程周期解的高階近似解析解，選用諧波平衡法是最簡(jiǎn)單有效的方法。將諧波平衡法求解非線性常微分方程的一般過程進(jìn)行了分解，將過程歸結(jié)為先設(shè)解為傅里葉級(jí)數(shù)形式并代入微分方程，經(jīng)過諧波平衡得到傅里葉系數(shù)方程組，最后求解方程組；

（2）對(duì)非線性項(xiàng)的傅里葉展開過程進(jìn)行了改進(jìn)，把多個(gè)運(yùn)算步驟歸結(jié)為兩次計(jì)算過程統(tǒng)一的矩陣運(yùn)算，使程序的編寫變得簡(jiǎn)潔，運(yùn)算比Ameer Hassan的方法更高效。

（3）非線性項(xiàng)的傅里葉展開計(jì)算結(jié)果與數(shù)值方法一致；誤差分析表明算法的精度主要依賴于諧波階數(shù)的選取。

[1]樓京俊.基于混沌理論的線譜控制技術(shù)研究［D］.武漢：海軍工程大學(xué)，2006.

[2]徐道臨，呂永建，周加喜.非線性隔振系統(tǒng)動(dòng)力學(xué)特性分析的FFT多諧波平衡法[J].振動(dòng)與沖擊，2012，31（22）：39-44.

[3]曹利，馮奇，張樂樂.雙非線性隔振器的雙層隔振系統(tǒng)模型的建立[J].噪聲與振動(dòng)控制，2008，28（2）：1-3.

[4]路純紅，白鴻柏，楊建春.超低頻非線性隔振系統(tǒng)的研究[J].噪聲與振動(dòng)控制，2010，30（4）：10-12.

[5]G.Duffing.Erzwungene Schwingungen bei ver?nderlicher EigenfrequenzundihretechnischeBedeutung[M].Braunschweig,Vieweg&Sohn,1918.

[6]劉延柱，陳立群.非線性振動(dòng)[M].北京：高等教育出版社，2001.

[7]劉式適，劉式達(dá).物理學(xué)中的非線性方程[M].北京：北京大學(xué)出版社，2000.

[8]王本利，張相盟，衛(wèi)洪濤.基于諧波平衡法的含Iwan模型干摩擦振子非線性振動(dòng)[J].航空動(dòng)力學(xué)報(bào)，2013，28（1）：1-9.

[9]呂永建.FFT多諧波平衡法及其應(yīng)用[D].長(zhǎng)沙：湖南大學(xué)，2012.

[10]Lau S L,Cheung Y K.Amplitude incremental variational principle for non-linear vibration of elastic systems[J].ASME,Journal of Applied Mechanics,1981:48:959-964.

[11]Mickens R E.An introduction to nonlinear oscillations [M].Cambridge University Press，Cambridge，1981.

[12]Hamdan M N,Burton T D.On the steady state response and stability of non-linear oscillators using harmonic balance[J].Journal of Sound and Vibration，1993，166 (2):255-266.

[13]Kawakami,Kobayashi.Computationofperiodic solutionsforpolynomialnonlinearequations[J].Electronics&Communications in Japan,1981,64-A(8): 7-14.

[14]Ameer Hassan.On the periodic and chaotic responses of duffing’s oscillator[D].Washington:Washington state university,1989.

[15]沈建和.非線性振動(dòng)系統(tǒng)的分岔、混沌及相關(guān)控制[D]．廣州：中山大學(xué)，2008.

[16]鈕耀斌，王中偉.基于橢圓函數(shù)諧波平衡法的二元機(jī)翼非線性顫振[J].工程力學(xué)，2013，30（4）：461-465.

[17]晏致濤，張海峰，李正良.基于增量諧波平衡法的覆冰輸電線舞動(dòng)分析[J].振動(dòng)工程學(xué)報(bào)，2012，25（2）：161-166.

[18]何漢林，梅家斌.數(shù)值分析[M].北京：科學(xué)出版社，2007.

[19]李厚儒，南敬昌.擬牛頓粒子群算法在非線性電路諧波平衡方程中的應(yīng)用[J].計(jì)算機(jī)應(yīng)用與軟件，2013，30（2）：103-105.

[20]李曉勇.基于強(qiáng)Duffing模型的隔振裝置混沌特性參數(shù)研究[D].廣州：廣州大學(xué)，2011.

[21]薛定宇，陳陽泉.高等應(yīng)用學(xué)習(xí)問題的MATLAB求解[M].北京：清華大學(xué)出版社，2004.

Simplification of Fourier Expansion in the High-order Harmonic Balance Method for Nonlinear Ordinary Differential Equations

TANG Yuan-zhang,WENG Xue-tao, LOU Jing-jun,LIN Xiong-wei,ZHANGHui

(1.College of Power Engineering,Naval University of Engineering,Wuhan 430033,China; 2.National Key Laboratory on Ship’s Vibration&Noise,Wuhan 430033,China)

A simplified computation method of the high-order harmonic solution for nonlinear ordinary differential equations is discussed.Fourier expansion procedure of the equation with quadratic or cubic terms is improved and simplified.The procedure consists of two steps of matrix operation with the same computation process so that the algorithm is easier to program than that of the previous equation.Results of the solution for the Duffing equation using this method show that the high-order harmonic solution is in good agreement with its numerical solution,but more efficient than the latter one.

vibration and wave;nonlinear ordinary differential equation;Duffing equation;Fourier expansion; harmonic balance method

O175.14

ADOI編碼：10.3969/j.issn.1006-1335.2014.02.007

1006-1355(2014)02-0028-06

2013-07-09

國(guó)家自然科學(xué)基金青年基金資助項(xiàng)目（基金編號(hào)：51009143）；

高等學(xué)校全國(guó)優(yōu)秀博士學(xué)位論文作者專項(xiàng)資金資助項(xiàng)目（基金編號(hào)：201057）

唐元璋（1985-），男，湖南湘潭人，碩士，從事非線性動(dòng)力學(xué)、振動(dòng)與噪聲控制研究。

E-mail:417006772@qq.com