，

(東北林業(yè)大學(xué) 理學(xué)院， 黑龍江 哈爾濱 150036)

雙曲型偏微分方程是描述振動(dòng)或波動(dòng)現(xiàn)象的一類(lèi)偏微分方程， 雙曲型偏微分方程的參數(shù)估計(jì)問(wèn)題廣泛存在于自然科學(xué)和工程技術(shù)領(lǐng)域中[1]， 因此研究雙曲型偏微分方程的參數(shù)估計(jì)方法在理論和實(shí)際方面都具有重要意義。 關(guān)于雙曲型方程的參數(shù)估計(jì)方法的研究已有一些進(jìn)展， 如差分演化算法[1-6]、 最佳攝動(dòng)量法[7-8]、 遺傳算法[9-10]等。 在理論上最完備而且行之有效的方法， 是由前蘇聯(lián)科學(xué)家Tikhonov以第一類(lèi)算子方程為基本框架于20世紀(jì)60年代初創(chuàng)造性提出且后來(lái)得到深入發(fā)展的正則化方法[11-12]， 但是它也有一定的局限性，例如正則化參數(shù)的選擇困難、 代步數(shù)的難以確定、 通用性的缺乏等。雖然雙曲線型方程的參數(shù)估計(jì)方法多種多樣，但是多數(shù)都處于理論計(jì)算階段，應(yīng)用到實(shí)際中的效果并不理想，并且一種方法大多都只能針對(duì)一類(lèi)問(wèn)題；因此研究新理論、 探索新方法是必要的，特別是通過(guò)采樣數(shù)據(jù)來(lái)估算方程中參數(shù)的方法頗具實(shí)際意義。本文中從新的角度出發(fā)，在僅已知模型類(lèi)型和觀測(cè)數(shù)據(jù)的條件下，將偏微分中的數(shù)值差分思想與最小二乘理論[13]結(jié)合，給出一種常系數(shù)非齊次雙曲型方程的參數(shù)估計(jì)方法，并通過(guò)實(shí)例模擬來(lái)驗(yàn)證該方法的可行性。

1 理論推導(dǎo)

含2個(gè)自變量x、y和未知函數(shù)u的二階線性常系數(shù)偏微分方程的一般形式為

Auxx+Buxy+Cuyy+Dux+Euy+Fu=g(x,y)，

(1)

式中：A、B、C、D、E、F為常數(shù)；g(x,y)為已知函數(shù)。當(dāng)B2-4AC>0時(shí)，該方程為雙曲型方程[14]。

將u(x,y)在區(qū)域[a,b]×[c,d](其中a、b、c、d均大于或等于0)分別取步長(zhǎng)h1和h2作2族與x軸、y軸平行的直線，得到矩形網(wǎng)格，2族直線的交點(diǎn)(a+ih1,c+jh2)稱為網(wǎng)點(diǎn)或節(jié)點(diǎn)[15]，記為(xi,yj)或(i,j)，節(jié)點(diǎn)處的函數(shù)值記為u(xi,yj)或uij(i,j=1,2,…)。

給定一組數(shù)據(jù)(u(k)(xi-1,yj),u(k)(xi,yj-1),u(k)(xi,yj),u(k)(xi,yj+1),u(k)(xi+1,yj),u(k)(xi+1,yj+1),g(k)(xi,yj)),k=1,2,…,N,N∈+,令

(2)

得到一組數(shù)據(jù)

(3)

方程(1)變?yōu)椴罘址匠?/p>

Fu(xi,yj)=g(xi,yj)。

(4)

方程(4)可以看作預(yù)測(cè)模型，在已知參數(shù)A、B、C、D、E、F的條件下，可以利用u(xi-1,yj)、u(xi,yj-1)、u(xi,yj)、u(xi,yj+1)、u(xi+1,yj)的值，估算u(xi+1,yj+1)的值。

將方程(1)看作多元線性函數(shù)，則把二階二維常系數(shù)雙曲型方程的參數(shù)估計(jì)問(wèn)題轉(zhuǎn)化成線性回歸模型的參數(shù)估計(jì)問(wèn)題，下面估計(jì)A、B、C、D、E、F的值。

設(shè)有N組數(shù)據(jù)(3)，利用最小二乘原理，令

Y=Xβ+e,

E(e)=0,

Y′Y=2Y′Xβ+β′X′Xβ,

Cov(e)=σ2In，

式中：σ為常數(shù)；In為n階單位矩陣，n∈+。

對(duì)β求偏導(dǎo)，并令其為0，得

X′Xβ=X′Y。

當(dāng)方程(1)的解u(x,y)為非二次及二次以下的代數(shù)多項(xiàng)式時(shí)，|X′X|一定不為0，因此X′X是可逆的，從而得到參數(shù)β的估計(jì)值

(5)

2 數(shù)值模擬與驗(yàn)證

為了驗(yàn)證上述方法的有效性，以雙曲型方程初值問(wèn)題

(6)

為例，對(duì)其參數(shù)進(jìn)行估計(jì)及數(shù)值模擬。

通過(guò)計(jì)算可知，上述雙曲型方程的解為

利用這N組數(shù)據(jù)，根據(jù)第1節(jié)中給出的方法來(lái)驗(yàn)證估算方程(6)中的參數(shù)，以此驗(yàn)證第1節(jié)中方法的有效性。

設(shè)有方程

Auxx+Buxy+Cuyy+Dux+Euy+Fu=xy，

(7)

利用N組數(shù)據(jù)估算A、B、C、D、E、F的估計(jì)值。

設(shè)h1=h2=h， 利用式(5)， 得出方程(7)中各參數(shù)的估計(jì)值， 見(jiàn)表1。 由表可知， 當(dāng)h→0時(shí)，A→1,B→3,C→2,D→0,E→0,F→0，因此方程(7)可以轉(zhuǎn)化為

Auxx+Buxy+Cuyy=xy。

(8)

表1 方程(7)中各參數(shù)的估計(jì)值

下面分3種情況討論方程(8)中各參數(shù)的估計(jì)值。

1)當(dāng)h1=h2=h時(shí)，利用式(5)，得到方程(8)中A、B、C的估計(jì)值， 見(jiàn)表2。 由表可知， 當(dāng)步長(zhǎng)h≤0.05時(shí)， 可以較準(zhǔn)確地估計(jì)出前面雙曲型方程(6)中的參數(shù)。 圖1所示為各參數(shù)估計(jì)值的相對(duì)誤差與步長(zhǎng)h的關(guān)系。 由圖可知，h越小， 方程各參數(shù)的相對(duì)誤差也越小， 當(dāng)h<0.05時(shí)，A的相對(duì)誤差小于0.4，B的相對(duì)誤差小于0.4，C的相對(duì)誤差小于0.2。這說(shuō)明第1節(jié)中給出的方法是有效的。

表2 在步長(zhǎng)h1=h2的條件下方程(8)中各參數(shù)的估計(jì)值

(a)參數(shù)A(b)參數(shù)B(c)參數(shù)C圖1 在步長(zhǎng)h1=h2的條件下各參數(shù)的相對(duì)誤差與步長(zhǎng)h的關(guān)系

2)固定h1=0.01，利用式(5)得出方程(8)中參數(shù)A、B、C的估計(jì)值，見(jiàn)表3。圖2所示為各參數(shù)的相對(duì)誤差與步長(zhǎng)h2的關(guān)系。由表3、圖2可知,當(dāng)固定h1=0.01時(shí)，并不是h2越小，A、B、C的相對(duì)誤差越小，而是h1和h2滿足一定關(guān)系時(shí)，A、B、C的估計(jì)值越接近真實(shí)值。當(dāng)h2是h1的2～3倍時(shí)，A、B、C的估計(jì)值更準(zhǔn)確。

表3 在步長(zhǎng)h1=0.01的條件下方程(8)中各參數(shù)的估計(jì)值

(a)參數(shù)A(b)參數(shù)B(c)參數(shù)C圖2 在步長(zhǎng)h1=0.01的條件下各參數(shù)的相對(duì)誤差與步長(zhǎng)h2的關(guān)系

3)固定h2=0.01，利用式(5)得出方程(8)中參數(shù)A、B、C的估計(jì)值，見(jiàn)表4。圖3所示為各參數(shù)的相對(duì)誤差與步長(zhǎng)h1的關(guān)系。由表4、圖3可知,當(dāng)固定h2=0.01時(shí)，h2越小,A、B、C的相對(duì)誤差也并非越小，而是h1和h2滿足一定關(guān)系時(shí)，A、B、C的估計(jì)值越接近真實(shí)值。

3 結(jié)論與討論

1)將數(shù)值差分思想與最小二乘理論結(jié)合，給出一種二階常系數(shù)非齊次雙曲型方程的參數(shù)估計(jì)方法，即利用依次采樣數(shù)據(jù)估算二階常系數(shù)非齊次雙曲型方程中的參數(shù)。如果雙曲型方程的解為非二次及低于二次的代數(shù)多項(xiàng)式，則本文中的方法具有可行性，并通過(guò)實(shí)例驗(yàn)證了該方法是可行的。

表4 在步長(zhǎng)h2=0.01的條件下方程(8)中各參數(shù)的估計(jì)值

(a)參數(shù)A(b)參數(shù)B(c)參數(shù)C圖3 在步長(zhǎng)h2=0.01的條件下各參數(shù)的相對(duì)誤差與步長(zhǎng)h1的關(guān)系

2)討論了不同步長(zhǎng)對(duì)參數(shù)的相對(duì)誤差的影響， 并得出結(jié)論： 如果步長(zhǎng)h1=h2=h， 則h越小， 其參數(shù)的相對(duì)誤差越?。?如果步長(zhǎng)h1≠h2， 則h1和h2需要滿足一定關(guān)系才能得出較好的結(jié)果。 由此可知， 步長(zhǎng)的選取與方程自身的特點(diǎn)密切相關(guān)。 如何根據(jù)方程自身特點(diǎn)選取有效的步長(zhǎng)是需要繼續(xù)探討的問(wèn)題。

3)文中的方法對(duì)于解為二次及低于二次的代數(shù)多項(xiàng)式并不適用，這是以后需要繼續(xù)探索的問(wèn)題。本文中以二階常系數(shù)非齊次雙曲型方程為研究對(duì)象，該方法也可以推廣到常系數(shù)非齊次的拋物型和橢圓型方程的參數(shù)估計(jì)中。如果進(jìn)一步討論，還可以估算出偏微分方程的邊界條件和初始條件，把這種方法推廣到實(shí)際應(yīng)用中。另外，本文中沒(méi)有考慮變系數(shù)非齊次偏微分方程的估計(jì)方法，如何準(zhǔn)確、有效地將該方法運(yùn)用到變系數(shù)偏微分方程的參數(shù)估計(jì)中是需要深入研究的問(wèn)題。

[1] 劉會(huì)超, 吳志健, 李煥哲, 等. 基于差分演化算法的雙曲型方程參數(shù)識(shí)別[J]. 武漢大學(xué)學(xué)報(bào)(理學(xué)版),2015,61(2): 117-123.

[2] 熊盛武, 李元香, 康立山， 等. 用演化算法求解拋物型方程擴(kuò)散系數(shù)的識(shí)別問(wèn)題[J]. 計(jì)算機(jī)學(xué)報(bào), 2000, 23(3): 261-265.

[3] 熊盛武, 李元香. 演化參數(shù)反演方法[J]. 武漢大學(xué)學(xué)報(bào)(理學(xué)版), 2001, 47(1): 37-41.

[4] STORN R, PRICE K. Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces[J]. Journal of Global Optimization, 1997, 11(4): 341-359.

[5] RAHNAMAYAN S, TIZHOOSH H R, SALAMA M M A. Opposition-based differential evolution[J]. IEEE Transactions on Evolutionary Computation, 2008, 12(1): 64-79.

[6] 吳志健. 演化優(yōu)化及其在微分方程反問(wèn)題中的應(yīng)用[D]. 武漢： 武漢大學(xué)， 2004.

[7] 彭亞綿, 楊?lèi)?ài)民, 龔佃選， 等. 改進(jìn)的最佳攝動(dòng)量法及在反問(wèn)題中的應(yīng)用[J]. 數(shù)學(xué)的實(shí)踐與認(rèn)識(shí), 2011, 41(5): 186-189.

[8] 彭亞綿. 雙曲型方程參數(shù)識(shí)別反問(wèn)題的解法[J]. 河北理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007, 29(3): 105-109.

[9] 張世梅. 二維偏微分方程反問(wèn)題的遺傳算法研究[D]. 西安： 西安理工大學(xué), 2005.

[10] 王小平, 曹立明. 遺傳算法[M]. 西安: 西安交通大學(xué)出版社, 2002: 22-68.

[11] RUDIN L I, OSHER S, FATEMI E. Nonlinear total variation based noise removal algorithms[J]. Physica: D: Nonlinear Phenomena, 1992, 60(1/2/3/4): 259-268.

[12] 吳建成, 張大力, 劉家琦. 一維波動(dòng)方程反問(wèn)題求解的正則化方法[J]. 計(jì)算物理, 1995, 12(3): 415-420.

[13] 王松桂, 陳敏, 陳立萍. 線性統(tǒng)計(jì)模型線性回歸與方差分析[M]. 北京: 高等教育出版社, 2014: 28-39.

[14] 謝鴻政, 楊楓林. 數(shù)學(xué)物理方程[M]. 北京: 科學(xué)出版社, 2008: 25-47.

[15] 李榮華, 劉播. 微分方程數(shù)值解法[M]. 北京: 高等教育出版社，2009: 67-76.