二維二階常系數(shù)雙曲型方程的參數(shù)估計方法

，

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

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

1 理論推導

含2個自變量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ù)。當B2-4AC>0時，該方程為雙曲型方程[14]。

將u(x,y)在區(qū)域[a,b]×[c,d](其中a、b、c、d均大于或等于0)分別取步長h1和h2作2族與x軸、y軸平行的直線，得到矩形網(wǎng)格，2族直線的交點(a+ih1,c+jh2)稱為網(wǎng)點或節(jié)點[15]，記為(xi,yj)或(i,j)，節(jié)點處的函數(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)可以看作預測模型，在已知參數(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ù)估計問題轉(zhuǎn)化成線性回歸模型的參數(shù)估計問題，下面估計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∈+。

對β求偏導，并令其為0，得

X′Xβ=X′Y。

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

(5)

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

為了驗證上述方法的有效性，以雙曲型方程初值問題

(6)

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

通過計算可知，上述雙曲型方程的解為

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

設(shè)有方程

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

(7)

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

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

Auxx+Buxy+Cuyy=xy。

(8)

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

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

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

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

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

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

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

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

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

3 結(jié)論與討論

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

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

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

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

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

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

[2] 熊盛武, 李元香, 康立山， 等. 用演化算法求解拋物型方程擴散系數(shù)的識別問題[J]. 計算機學報, 2000, 23(3): 261-265.

[3] 熊盛武, 李元香. 演化參數(shù)反演方法[J]. 武漢大學學報(理學版), 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)化及其在微分方程反問題中的應用[D]. 武漢： 武漢大學， 2004.

[7] 彭亞綿, 楊愛民, 龔佃選， 等. 改進的最佳攝動量法及在反問題中的應用[J]. 數(shù)學的實踐與認識, 2011, 41(5): 186-189.

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

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

[10] 王小平, 曹立明. 遺傳算法[M]. 西安: 西安交通大學出版社, 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] 吳建成, 張大力, 劉家琦. 一維波動方程反問題求解的正則化方法[J]. 計算物理, 1995, 12(3): 415-420.

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

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

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