陳風(fēng)雷，馬正義,2

(1.浙江理工大學(xué)理學(xué)院，杭州 310018；2.麗水學(xué)院工程與設(shè)計(jì)學(xué)院，浙江麗水 323000)

?

熱方程邊界值決定反問題的數(shù)值方法

陳風(fēng)雷1，馬正義1,2

(1.浙江理工大學(xué)理學(xué)院，杭州 310018；2.麗水學(xué)院工程與設(shè)計(jì)學(xué)院，浙江麗水 323000)

給出了求解熱方程邊界值決定反問題的直接差分法和配置法兩種數(shù)值算法。直接差分法采用Euler向后差分格式；配置法把反問題化為擬解問題，利用Lagrange插值基函數(shù)構(gòu)造有限維逼近，從而將問題轉(zhuǎn)化為求解一個(gè)代數(shù)方程組，最后采用正則化方法求解。數(shù)值結(jié)果表明：直接差分法的反演結(jié)果在區(qū)間的左端有較強(qiáng)的振蕩，而配置法能夠整體恢復(fù)左邊界值，并且配置法可靈活選擇正問題的數(shù)值格式獲得較高的精度。

邊界值決定反問題；直接差分法；配置法；正則化； Lagrange插值

0 引 言

自然科學(xué)與工程技術(shù)上的許多問題都可以歸結(jié)為反問題，其中一類重要的反問題稱為邊界值決定反問題，即由可接觸的邊界或容易接觸的邊界(值)的信息決定不可接觸或不易接觸的邊界(值)的信息。例如，織物外側(cè)的溫度和濕度容易測量，但織物與人體之間的溫度不易測量。類似的例子還有無損檢測、高爐煉鐵等，甚至地震預(yù)測也可以歸為這類問題[1-2]。

本文研究熱方程的邊界值決定反問題。該問題是不適定的，輸入數(shù)據(jù)的微小誤差可能導(dǎo)致數(shù)值解的巨大變化，因此需要采用正則化技巧給出穩(wěn)定的數(shù)值方法[3]。目前已有不少這方面的工作，主要可以歸為兩類方法：差分法和Fourier變換方法。

對于差分法，Murio[3]提出了完全顯示且穩(wěn)定的空間差分推進(jìn)格式。利用具有平均核的離散卷積過濾噪聲數(shù)據(jù)，Murio把該類反問題轉(zhuǎn)化為一個(gè)適定性問題，然后使用空間推進(jìn)的差分格式求解。Elden[4]研究了時(shí)間差分推進(jìn)格式的正則化效應(yīng)，發(fā)現(xiàn)時(shí)間差分可以防止解在高頻部分發(fā)生爆破，即時(shí)間方向的差分具有正則化效應(yīng)。

另外一類為Fourier變換方法，該方法把偏微分方程從空間域轉(zhuǎn)化到頻率域，從而獲得復(fù)變量的常微分方程。對該常微分方程，Hao[5]給出了Fourier逆變換下的顯示解，發(fā)現(xiàn)反問題不穩(wěn)定正是高頻部分造成的。對顯示解的不同處理可獲得不同的數(shù)值方法，如Fourier正則化方法、譜正則化方法、Tikhonov法、小波方法和最優(yōu)過濾法等[6-9]，這些方法統(tǒng)稱為Fourier變換方法。

以上兩類方法各有優(yōu)缺點(diǎn)。差分方法簡單，容易實(shí)施，但對復(fù)雜方程很難給出誤差估計(jì)，而且不同差分格式的結(jié)果往往呈現(xiàn)較大差異。Fourier變換方法對空間變量進(jìn)行Fourier變換，這要求方程系數(shù)僅依賴于空間變量，尤其對非線性問題，F(xiàn)ourier變換難以實(shí)施，這限制了Fourier變換方法的適用范圍。

Hasanov等[10]用配置法求解了一類反問題，即時(shí)間后向的拋物型問題(backward parabolic problems，BPP)，獲得了穩(wěn)定且較高精度的數(shù)值解。本文將配置法用于求解熱方程邊界值決定反問題。配置法是對原問題的有限維逼近，通過對反演函數(shù)進(jìn)行插值，將反問題轉(zhuǎn)化為一系列以插值基函數(shù)為輸入數(shù)據(jù)的正問題和一個(gè)不適定的代數(shù)方程組，而正問題是穩(wěn)定的，因此配置法可以有效地抑制反問題的不穩(wěn)定性。本文給出了直接差分法和配置法的數(shù)值結(jié)果，并對兩種方法進(jìn)行了比較。

1 熱方程邊界值決定反問題

熱方程邊界值決定反問題可用下面一般的方程表示：

(1)

其中：u∈C2,1((0,1)×(0,T))為溫度；κ(x,t)∈C1,0((0,1)×(0,T))表示熱傳導(dǎo)率；Θ(x,t)表示源項(xiàng)(輻射項(xiàng)和蒸發(fā)冷凝項(xiàng))；v(x,t)為初始條件；φ(t)，Q(t)表示該方程邊界條件； T表示時(shí)間變量的上界。反問題的目的是通過求解問題(1)來決定左邊界值u(0,t)。本文稱問題(1)為IP(inverse problems)。

上述反問題可等價(jià)描述為求解P(t)，使得滿足方程組：

(2)

且

(3)

易知，對于可測量數(shù)據(jù)：初邊值、熱傳導(dǎo)率κ(x,t)和源項(xiàng)Θ(x,t)，若初邊值發(fā)生變化時(shí)，若問題(2)的解并不是線性依賴于左邊界值P(t)的解，但可以通過如下處理進(jìn)行線性化。令u=w+z，其中w、z分別滿足：

顯然w關(guān)于P(t)是線性的，而z可以直接求解。

則反問題歸結(jié)為求算子方程

f(P)=0

的零點(diǎn)?？煽紤]最小二乘解，即求如下泛函的極小元：

F:L2[0,T]→[0,+∞)，

(4)

以上優(yōu)化問題的解稱為IP的擬解。

2 數(shù)值算法

本節(jié)考慮IP的數(shù)值算法，即直接差分法和配置法。

2.1 直接差分法

(5)

令s=τ/h2，式(5)可等價(jià)如下線性方程組：

即：

(6)

其中：

2.2 配置法

用配置法求解(4)的算法如下：

Step 1 用Lagrange基函數(shù)對邊界值P(t)進(jìn)行插值，即：

由于所考慮的問題是線性的，因此有：

Step 2 泛函(4)可寫為

由極值的必要性：

從而有線性方程組

Aλ=b

(7)

其中A=(aij)，b=(b0,b1,…,bn)T，而

Step 4 如果矩陣A的條件數(shù)較大，則采用正則化方法求解，如廣義交叉校驗(yàn)Tikhonov正則化法、廣義交叉校驗(yàn)TSVD正則化法[11-13]，否則直接用最速下降法[14]求解。

3 數(shù)值模擬

算例1:考慮如下問題：

易知該問題有唯一解u(x,t)=ex+t，從而u(0,t)=et。

a)直接差分法的數(shù)值結(jié)果

取步長τ=1/100，h=1/200，直接差分法的數(shù)值解與精確解見圖1和圖2。

從圖1和圖2中可以得出，左邊界值在區(qū)間的左側(cè)具有明顯的振蕩，為此進(jìn)行截?cái)?，見圖3。

(a)精確解

(b)數(shù)值解圖1 u(x,t)精確解與數(shù)值解

圖2 u(0,t)數(shù)值解與精確解

圖3 u(0,t)截?cái)嗪蟮臄?shù)值解與精確解

圖3表明，截?cái)嗪蟮膱D像與精確解吻合。在實(shí)際運(yùn)用中，若考慮長時(shí)間后的擴(kuò)散或傳遞行為，差分法不論在計(jì)算時(shí)間上還是數(shù)值精度上都滿足工業(yè)要求。

b)配置法的數(shù)值結(jié)果

取Lagrange基函數(shù)的階n=5，反問題離散為N=100,M=200，計(jì)算(7)時(shí)矩陣A的條件數(shù)為1.594×104，可直接采用最速下降法求解。

對右邊界u(L,T)數(shù)據(jù)在不同噪聲水平下，u(0,t)數(shù)值解與精確解如圖4所示。

圖4 不同噪聲δ下u(0,t)數(shù)值解與精確解

圖4表明：配置法的反演結(jié)果并沒有產(chǎn)生振蕩，且具有較高的精度，因此對整體反演結(jié)果，配置法比較合適。

表1 配置法誤差結(jié)果

通過不同階的Lagrange基函數(shù)進(jìn)行數(shù)值計(jì)算，結(jié)果如圖5所示。

圖5 噪聲δ=1%時(shí)u(0,t)的數(shù)值解與精確解

從圖5(a)可以發(fā)現(xiàn)數(shù)值解誤差比較小，但在提高Lagrange基函數(shù)的階數(shù)為15時(shí)(見圖5(b))，計(jì)算精度并沒有因?yàn)殡A的增大而提高，這可能是高次插值的龍格現(xiàn)象[14]。為此可嘗試采用Chebyshev多項(xiàng)式的零點(diǎn)作為插值節(jié)點(diǎn)以避免龍格現(xiàn)象，數(shù)值結(jié)果如圖6所示。從圖6中可以發(fā)現(xiàn)，在端點(diǎn)處同樣出現(xiàn)較大誤差，這是算法本身的原因。事實(shí)上，對以每個(gè)基函數(shù)為左邊界值的正問題，誤差在原點(diǎn)附近都是最大，因此誤差是由正問題的數(shù)值算法導(dǎo)致的。

圖6 噪聲δ=1%時(shí)u(0,t)的數(shù)值解和精確解

算例2 考慮如下方程：

已知上述方程精確解：u(x,t)=exsint，從而初始條件為：u(0,t)=sint.

a)差分法的數(shù)值結(jié)果

取步長為τ=1/100，h=1/400邊界u(0,t)數(shù)值解與精確解如圖7所示。

圖7中同樣發(fā)現(xiàn)：在若干個(gè)時(shí)刻的數(shù)值結(jié)果有明顯的振蕩，并且在振蕩后數(shù)值解與真實(shí)解非常接近。

b)配置法的數(shù)值結(jié)果

取Lagrange基函數(shù)的次數(shù)n=5，反問題空間離散取N=100，M=200，計(jì)算(7)時(shí)其矩陣的條件數(shù)為86.3437，直接采用最速下降法。

圖7 u(0,t)數(shù)值解與精確解

圖8 噪聲δ=0.01時(shí)u(0,t)數(shù)值解與真實(shí)解

從圖8可以看出，配置法克服了直接差分法在起始時(shí)刻產(chǎn)生的較強(qiáng)的振蕩，從而獲得整個(gè)區(qū)間上的穩(wěn)定解。

4 結(jié) 論

本文采用直接差分法和配置法研究了一類熱方程邊界值決定反問題，并給出相關(guān)數(shù)值計(jì)算。直接差分法只需簡單迭代，用時(shí)短，但數(shù)值結(jié)果表明，直接差分法在區(qū)間左端具有明顯的振蕩。因此，若考慮長時(shí)間后的擴(kuò)散或傳遞行為，直接差分法在時(shí)間損耗和數(shù)值精度上都能滿足工業(yè)要求。配置法把反問題轉(zhuǎn)化為一系列的正問題，通過靈活選取正問題的數(shù)值格式，不僅能夠獲得較高的精度，而且對于初邊值在不同的噪聲水平下都可以避免較強(qiáng)的振蕩。但是配置法需要求解大量正問題，為此可采用并行計(jì)算縮短計(jì)算時(shí)間。綜上分析，配置法對熱方程邊界值決定反問題比直接差分法具有更強(qiáng)的適用性。

[1] 徐定華.紡織材料熱濕傳遞數(shù)學(xué)模型及設(shè)計(jì)反問題[M].北京:科學(xué)出版社,2014:1-12.

[2] 戴會(huì)超,王玲玲.工程水力學(xué)反問題研究及應(yīng)用[J].四川大學(xué)學(xué)報(bào),2006,38(1):15-19.

[3] MURIO D A. The mollification method and the numerical solution of the inverse heat conduction problem by finite differences [J]. Computers & Mathematics with Applications,1989,17(10):1385-1396.

[4] ELDEN L. Numerical solution of the sideways heat equation bydifference approximation in time [J]. Inverse Problem,1995,11(4):913-923.

[5]HAO D N. On a sideways parabolic equation[J]. Inverse Problems,1998,13(2):297-309.

[6] CARASSO A. Determining surface temperatures from interior observations[J]. SIAM JAppl Math,1982,42(3):558-574.

[7] CHOUDHURYAH. Wavelet Method for Numerical Solution of Parabolic Equations[J/OL]. Journal of Computational Engineering,2014:1-12.http://dx.doi.org/10.1155/2014/346731.

[8] SEIDMAN T I,ELDEN L. An optimal filtering method for the sideways heat equation[J].Inverse Problems,1990,6(4):681-696.

[9] TAUTENHAHN U. Optimal stable approximations for the sideways heat equation[J]. Inverse Ill-Posed Problems,1997,5(3):287-307.

[10] HASANOVA,MUELLER J L .A numerical method for backward parabolic problems with non-self adjoint elliptic operators[J].Applied Numerical Mathematics,2001,37(1/2):55-78.

[11] 劉繼軍.不適定問題的正則化方法及應(yīng)用[M].北京:科學(xué)出版社,2005:46-55.

[12] 韓波,李莉.非線性不適定問題的求解方法及應(yīng)用[M].北京:科學(xué)出版社,2011:25-118.

[13] HANSEN PC. The truncated SVD as a method for regularization [J]. Bit Numerical Mathematics,2010,27(4):534-553.

[14] 李慶揚(yáng),王能超,易大義.數(shù)值分析[M].北京:清華大學(xué)出版社,2012:202-204.

(責(zé)任編輯： 康 鋒)

Numerical Methods of Inverse Problem of Boundary Value Determination of Heat Equation

CHENFenglei1,MAZhengyi1,2

(1.School of Science, Zhejiang Sci-Tech University, Hangzhou 310018, China； 2.Institute of Engineering and Design, Zhejiang Lishui University, Lishui 323000, China)

In this paper, a direct difference method and a collocation method are applied to solve an inverse problem of boundary value determination of heat equation. The direct difference method adopts the Euler backward difference scheme. The collocation method first formulates the inverse problem to a quasi-solution problem and then uses base function of Lagrange interpolation to construct finite dimension approximation. In this way, the original problem is turned into solution of a system of algebraic equations. Finally, regularization method is used to solve. Numerical result shows that the inversion result of direct difference method presents strong fluctuation at the left end of interval, while the collocation method can recover the left boundary value on the whole and is easy to obtain higher precision by choosing appropriate numerical format of the direct problem.

inverse problem of boundary value determination; direct difference method; collocation method; regularization; Lagrange interpolation

10.3969/j.issn.1673-3851.2016.11.025

2016-04-16

國家自然科學(xué)基金項(xiàng)目(11071221)；浙江省自然科學(xué)基金項(xiàng)目(LY14A010005)

陳風(fēng)雷(1990-)，男，安徽淮北人，碩士研究生，主要從事數(shù)理方程的反問題方面的研究。

馬正義，E-mail：ma-zhengyi@163.com

O241.5； O241.82

A

1673- 3851 (2016) 06- 0945- 06 引用頁碼： 110802