一種求解非良態(tài)高階ODE初值問題的類解析法

陳 毅,伍儒康,王 偉,凌家勝

(南京工程學(xué)院機(jī)械工程學(xué)院, 江蘇 南京 211167)

1 常微分方程的研究現(xiàn)狀

微分方程的解是存在且唯一的,但由于理論方法的局限性,很多高階微分方程無法求出解析解.從實(shí)際應(yīng)用上講,人們往往不需要解析解,而是關(guān)心某個(gè)定義范圍內(nèi)求出的對(duì)應(yīng)近似值,即數(shù)值解[1-2].

常微分方程(簡(jiǎn)稱ODE)的數(shù)值解法包括歐拉方法、龍格-庫(kù)塔方法和差分法等,其中高階(包括三階)常微分方程一般可以化為一階方程組的初值問題進(jìn)行求解或者采用狀態(tài)空間法運(yùn)用經(jīng)典的4級(jí)四階龍格-庫(kù)塔方法進(jìn)行求解[3-5].但是差分法和龍格-庫(kù)塔方法都存在求解的局限性,當(dāng)微分方程三階以上時(shí),差分后的離散方程難以保證對(duì)角占優(yōu)、絕對(duì)收斂;使用龍格-庫(kù)塔方法求解非良態(tài)高階常微分方程時(shí),存在計(jì)算周期較長(zhǎng)的問題[6-8].本文針對(duì)上述問題,結(jié)合微分方程解析解法與數(shù)值解法,研究微分方程的類解析法.

2 三階ODE初值問題

定義三階ODE的初值問題:

(1)

式中,y為輸出,需要求解.

旋轉(zhuǎn)機(jī)械發(fā)電裝置導(dǎo)出的數(shù)學(xué)模型為:

(2)

式中:g(x)為輸入;p、q、s′和d分別為三階、二階、一階和y的系數(shù);p≠0且g(x)≠0,即對(duì)應(yīng)的三階微分方程是非齊次的.

3 類解析法的介紹

3.1 輸入函數(shù)的分解

周期為T的輸入函數(shù)g(x)按三角形式的傅里葉級(jí)數(shù)展開為:

(3)

式中,a0、ai、bi為傅里葉系數(shù).

解方程可得:

py?+qy″+s′y′+dy=

(4)

展開有:

(5)

解出每一個(gè)方程對(duì)應(yīng)的解析解,然后將上述2n+1項(xiàng)線性疊加.

3.2 特征方程的求解

式(4)對(duì)應(yīng)的特征方程可以通過Matlab的solve命令求得,也可以通過迭代得到.特征根的情況分為四種:

情況1,px3+qx2+s′x+d=p(x+r1)(x+r2)×(x+r3),其中r1、r2、r3為特征方程的3個(gè)全不相等實(shí)根;

情況2,px3+qx2+s′x+d=p(x+r1)[(x+α)2+β2],1個(gè)實(shí)根、2個(gè)共軛復(fù)根;

情況3,3個(gè)相等實(shí)根;

情況4,3個(gè)實(shí)根,包含2個(gè)相等實(shí)根.

本文主要討論情況1和情況2,情況3和情況4與情況1求解類似,且求解更簡(jiǎn)單,不再贅述.

3.3 拉氏變換與拉氏反變換

3.3.1 情況1

1) 輸入為e1=a0/2,利用拉氏變換中微分定理[9]可得:

py?+qy″+s′y′+dy=e1

(6)

變換為:

(7)

式中:e1為傅里葉多項(xiàng)式中的常數(shù)項(xiàng)a0/2;s為復(fù)數(shù)σ+jω;s′為微分方程一階系數(shù),即特征方程一次系數(shù);r1、r2、r3為情況1的3個(gè)全不相等實(shí)根;k1、k2、k3、k4為分母連乘真分式化成部分真分式的系數(shù).

將式(7)進(jìn)行通分得到s的三次方程,s的三次、二次和一次系數(shù)應(yīng)為0,常數(shù)項(xiàng)應(yīng)等于1,4個(gè)方程可確定4個(gè)系數(shù):

(8)

下面系數(shù)求解同理,不再贅述.通過拉氏反變換得到常數(shù)項(xiàng)輸入對(duì)應(yīng)的輸出為:

(9)

(10)

系數(shù)可以通過Matlab矩陣求解得到:

(11)

通過拉氏反變換可得:

k3e-r2t+k4e-r3t)

(12)

(13)

式中:ei+n+1為第i個(gè)余弦輸入的振幅,即bi;k11、k12、k2、k3、k4為部分真分式的系數(shù).

系數(shù)可以通過Matlab矩陣求解得到:

(14)

通過拉氏反變換可得:

k3e-r2t+k4e-r3t)

(15)

式(15)與式(12)的區(qū)別在于求解的系數(shù)不同.

3.3.2 情況2

(16)

系數(shù)的求解為:

(17)

對(duì)應(yīng)的拉氏反變換為:

(18)

(19)

系數(shù)通過Matlab矩陣求解得到：

(20)

同理,通過拉氏反變換可得:

(21)

3) 輸入為ei+n+1coswit(i=1,2,…,n),其中

(22)

系數(shù)通過Matlab矩陣求解得到：

(23)

拉氏反變換為:

(24)

3.4 最終輸出

將2n+1項(xiàng)線性疊加,得到最終的輸出函數(shù)方程為:

(25)

當(dāng)n較小時(shí),這是一個(gè)有意義的解析解;當(dāng)n較大時(shí),可以根據(jù)步長(zhǎng)得出數(shù)值解.

4 類解析法的運(yùn)用

求解

式中:右邊為系統(tǒng)輸入;左邊時(shí)變函數(shù)y(t)為系統(tǒng)輸出.

已知系統(tǒng)輸入為周期函數(shù)且周期為2.5 s,如圖1所示,圖1中僅展示了系統(tǒng)輸入的5個(gè)周期.

圖1 系統(tǒng)輸入隨時(shí)間變化圖

很明顯,上述三階ODE各系數(shù)數(shù)量級(jí)相差較大,采用本文類解析法進(jìn)行求解，步驟為:

1) 將系統(tǒng)輸入進(jìn)行傅里葉級(jí)數(shù)展開,展開階數(shù)n=4,從圖1中不難看出,9項(xiàng)傅里葉展開式可以很好地代替輸入函數(shù);

2) 運(yùn)用Matlab進(jìn)行求解特征根,求解結(jié)果為1個(gè)實(shí)根和2個(gè)共軛虛根;

3) 利用循環(huán)指令和矩陣運(yùn)算求出常數(shù)輸入的4個(gè)系數(shù)、正弦輸入的20個(gè)系數(shù)和余弦輸入的20個(gè)系數(shù),共計(jì)44個(gè)常數(shù),再回代到每個(gè)輸出分量的表達(dá)式中,得出每項(xiàng)輸入對(duì)應(yīng)的輸出;

4) 將輸出分量線性疊加,得出系統(tǒng)輸出,如圖2所示,可知本方法收斂;

圖2 系統(tǒng)輸出隨時(shí)間變化圖

5) 利用得到的輸出數(shù)據(jù)進(jìn)行一階差分、二階差分和三階差分,代入例題中與圖1的系統(tǒng)輸入進(jìn)行比較,明顯可見誤差較小,驗(yàn)證了本方法求解結(jié)果的正確性.

5 結(jié)語(yǔ)

本文對(duì)三階常系數(shù)ODE的數(shù)值求解方法進(jìn)行探討,給出了一個(gè)行之有效的求解方法——類解析法,即將三階常系數(shù)ODE的輸入按傅里葉級(jí)數(shù)展開成若干項(xiàng),將一個(gè)輸入的微分方程分解為若干輸入分量的微分方程,分別求解并線性疊加,減少了單個(gè)方程的求解難度.通過計(jì)算可以證明,該方法適應(yīng)性強(qiáng)、構(gòu)造簡(jiǎn)單、可復(fù)制性強(qiáng)、求解準(zhǔn)確且計(jì)算周期短,可以應(yīng)用于機(jī)械、液壓和電氣等三階系統(tǒng)中,是一個(gè)有效求解類解析解的數(shù)值解法.