李麥俠, 曲小鋼

(西安建筑科技大學(xué) 理學(xué)院， 陜西 西安 710055)

0 引言

本文采用直線法(Method of Lines,簡(jiǎn)稱MOL)對(duì)空間域進(jìn)行離散，并利用四階龍格-庫(kù)塔方法可以長(zhǎng)時(shí)間精確計(jì)算的特點(diǎn)，實(shí)現(xiàn)了對(duì)梁振動(dòng)問(wèn)題高精度、穩(wěn)定地求解.

1 MOL數(shù)值解法

設(shè)在矩形區(qū)域G{0≤x≤1;0≤t≤T}內(nèi)研究如下梁振動(dòng)方程的初邊值問(wèn)題.

其中,φ1(x)、φ2(x)、g1(t)、g2(t)、g3(t)、g4(t)

是給定函數(shù)[4]，φ1(x)、φ2(x)為初值,g1(t)、g2(t)、g3(t)、g4(t)為邊界值.

采用直線法對(duì)空間導(dǎo)數(shù)進(jìn)行離散[5,6].

圖1 空間變量x的剖分圖

假定問(wèn)題(Ⅰ)有充分光滑解u(x,t),可按直線法思路來(lái)求解.在方程(Ⅰ)中令x=xk,并利用如下二階中心差商：

2u(xk+1,t)+4u(xk,t)-2u(xk-1,t)+u(xk,t)-2u(xk-1,t)+u(xk-2,t)]=

4u(xk-1,t)+u(xk-2,t)]

再令uk(t)=u(xk,t),可得方程：

o(Δx4) (k=2,…,N-2)

(5)

(6)

利用初始條件(2)可得：

(7)

利用邊界條件(3)可得：

U0(t)=U(0,t)=g1(t)UN(t)=U(1,t)=g2(t)

(8)

利用邊界條件(4)及一階向前差商可得：

g3(t)

(9)

g4(t)

(10)

具有邊界條件 (8)、(9)、(10)的方程組(6)是以o(Δx4)的精確度逼近問(wèn)題(Ⅰ),此方程組稱為直線法方程組.(7)、(8)、(9)、(10)式共同構(gòu)成了常微分方程組的初邊值問(wèn)題,要得到此方程組的解析解u(xk,t)(k=1,2,…,N-1)相當(dāng)困難.

由于常微分方程初邊值問(wèn)題數(shù)值解的算法已相當(dāng)成熟,其中龍格-庫(kù)塔法[7]被證明是具有精度高、收斂、穩(wěn)定(在一定條件下)、計(jì)算過(guò)程中可以改變步長(zhǎng)、不需要計(jì)算高階導(dǎo)數(shù)等優(yōu)點(diǎn)的一種方法,故本文采用四階龍格-庫(kù)塔法求解以上的常微分方程邊值問(wèn)題的近似解u(xk,tj).

下面考慮更一般的一階常微分方程組的初邊值問(wèn)題：

其四階龍格-庫(kù)塔法數(shù)值求解公式[8]為：

yi,m+1=yim+(h/6)(ki1+2ki2+2ki3+ki4)

其中

(i=1,2,…,n;m=0,1,2,…)

其中，yim是第i個(gè)函數(shù)yi在節(jié)點(diǎn)tm=t0+mh處的近似值,h為積分步長(zhǎng).

根據(jù)上面的計(jì)算過(guò)程及公式，就可以進(jìn)行計(jì)算機(jī)編程，其具體步驟如下：

(1)對(duì)直線x=xk時(shí)的方程，用二階中心差商代替x的偏導(dǎo)數(shù)，得到線性常微分方程組；

(2)Matlab編程,用龍格-庫(kù)塔方法求解常微分方程組，得到數(shù)值解；

(3)調(diào)出數(shù)值解且與真值比較[9]，得到誤差分析.

按以上的步驟，便可求出線性梁振動(dòng)方程的數(shù)值解.

2 上機(jī)實(shí)現(xiàn)

利用Matlab語(yǔ)言編寫了上述線性二階常微分方程組邊值問(wèn)題數(shù)值解的程序,從而可得到u在t所屬區(qū)間[0≤t≤T]上的節(jié)點(diǎn){t0=0,t1=Δt,…,tj=jΔt,…,tM=MΔt}上的近似解u(xk,tj)，其中Δt=T/M.

3 收斂性及誤差分析

設(shè)rk(t)=u(xk,t)-Uk(t)，由Taylor展開(kāi)式[10,11]得到rk(t)的方程組

(11)

(12)

r2(t)-r1(t)+r0(t)=0,

rN(t)-2rN-1(t)+rN-2(t)=0

(13)

其中,k=0,1,2,…,N;此時(shí)0≤x≤1,0≤t≤T.

假定下式[12,13]

(14)

由(11)、(12)、(13)和Cauchy不等式[6]、微分(14) ，得到：

因?yàn)镮(0)=0，所以

(15)

其次，根據(jù)r0(t)=rN(t)=0，得到：

因此，由Cauchy不等式

最終得到如下估計(jì)：

4 結(jié)束語(yǔ)

本文給出了求梁振動(dòng)問(wèn)題的一種新的方法，即直線法.它具有精度高,編程簡(jiǎn)單,易于計(jì)算機(jī)實(shí)現(xiàn)等特點(diǎn),特別適用于非線性方程的數(shù)值求解,是一種實(shí)用有效的數(shù)值求解方法.

[1] 倪 平，高儀新.解梁的振動(dòng)方程的廣義方法(I)[J].東北師范大學(xué)學(xué)報(bào)(自然科學(xué)版)，1995(4):14-19.

[2] 曾文平.四階桿振動(dòng)方程的含參數(shù)四層顯式差分格式[J].華僑大學(xué)學(xué)報(bào)(自然科學(xué)版)，2002,23(2):116-121.

[3] 曾文平.解四階桿振動(dòng)方程新的兩類隱式差分格式[J].華僑大學(xué)學(xué)報(bào)(自然科學(xué)版)，2003,24(2):136-142.

[4] 鄒 佩，曲小鋼.梁振動(dòng)問(wèn)題的擬小波-精細(xì)時(shí)程積分法[J].陜西科技大學(xué)學(xué)報(bào)(自然科學(xué)版)，2011,29(6):140-144.

[5] 彭亞綿，閔 濤，張世梅，等.Burgers方程的MOL數(shù)值解法[J].西安理工大學(xué)學(xué)報(bào)，2004,20(3):276-280.

[6] 南京大學(xué)數(shù)學(xué)系.偏微分方程數(shù)值解法[M].北京:科學(xué)出版社，1979：393-417.

[7] 易晉生，顧安邦，王小松.基于MATLAB的公路橋梁車橋耦合數(shù)值計(jì)算方法[J].重慶交通大學(xué)學(xué)報(bào)(自然科學(xué)版)，2012,31(6):1 101-1 104.

[8] 毋玉芝.四階龍格-庫(kù)塔算法的C語(yǔ)言實(shí)現(xiàn)[J].焦作大學(xué)學(xué)報(bào),2001,15(1):55-57.

[9] 單雙榮.梁振動(dòng)方程的多辛Fourier擬譜算法[J].華僑大學(xué)學(xué)報(bào)(自然科學(xué)報(bào))，2006,27(3):234-237.

[10] 許士菊，王長(zhǎng)華.梁振動(dòng)方程的一個(gè)穩(wěn)定的有限差分近似[J].吉林化工學(xué)院學(xué)報(bào)，2007,24(1):79-81.

[11] 周良強(qiáng)，陳予怒，陳芳啟.梁振動(dòng)方程的多參數(shù)高精度格式[J].安慶師范學(xué)院學(xué)報(bào)(自然科學(xué)版)，2010,16(2):62-65.

[12] 戴嘉尊，邱建賢.微分方程數(shù)值解法[M].南京:東南大學(xué)出版社，2002.

[13] 孫志忠.偏微分方程數(shù)值解法[M].北京:科學(xué)出版社，2005.