唐錦萍

(黑龍江大學(xué) 數(shù)據(jù)科學(xué)與技術(shù)學(xué)院，哈爾濱150080)

1 引 言

在實(shí)際問題中，有很多數(shù)學(xué)物理問題的求解，都可以歸結(jié)為如下形式的第一類Fredholm積分方程的求解問題：

(1)

其中，A(x，y)，g(y)假定為連續(xù)有界函數(shù)，f(x)為待求解的未知函數(shù).方程(1)的求解，積分方程的離散是關(guān)鍵.根據(jù)離散方式的不同，積分方程的離散策略可以分為兩類，一類是采用數(shù)值積分公式，如復(fù)化梯形公式等，直接在離散節(jié)點(diǎn)上對(duì)方程(1)進(jìn)行離散[1]；另一類是先對(duì)區(qū)間分段，并給出未知函數(shù)f(x)的逼近函數(shù)，如分段Hermit插值函數(shù)，由于每個(gè)區(qū)間上逼近函數(shù)的解析形式是給定的，因此(1)式左端的積分在每一個(gè)小區(qū)間上可以被解析的計(jì)算出來，進(jìn)而將所有區(qū)間上的積分整合到一起，就可以得到方程(1)的離散形式[2].本文將采用第二類離散方式，利用三次樣條插值函數(shù)近似f(x).三次樣條函數(shù)具有分段、低階、光順的特性，經(jīng)常被用來作為光滑函數(shù)的近似[3-4].

由于第一類Fredholm積分方程的求解通常是不適定的，因此常采用正則化策略來克服，進(jìn)而得到穩(wěn)定的數(shù)值解.注意到在以往的利用三次樣條函數(shù)進(jìn)行積分方程離散的文獻(xiàn)中[5]，作者一般采用的是先正則化后離散的方式.這種方式先將積分方程的求解轉(zhuǎn)化為帶約束的優(yōu)化問題：約束條件就是積分方程(1)，優(yōu)化目標(biāo)是一個(gè)極小化泛函，通?？梢詾槲粗瘮?shù)f(x)的L2范數(shù)，或者其一階導(dǎo)數(shù)f′(x)、二階導(dǎo)數(shù)f″(x)的L2范數(shù)，其中后兩種優(yōu)化目標(biāo)能起到光滑化的作用.然后將帶約束的優(yōu)化問題轉(zhuǎn)化成無約束優(yōu)化問題，并用三次樣條函數(shù)離散，這樣就得到了關(guān)于未知函數(shù)在離散節(jié)點(diǎn)上的函數(shù)值的線性方程組.該線性方程組的解通常是存在的，且其解是穩(wěn)定的.與先正則化后離散的策略不同，本文將采用先離散后正則化的思路進(jìn)行求解，先對(duì)f(x)利用三次樣條函數(shù)進(jìn)行插值逼近，得到方程(1)的離散化線性方程組的形式后[6]，通過引入未知解的離散化光滑條件，將方程(1)的求解轉(zhuǎn)化為求解關(guān)于線性方程組的帶約束優(yōu)化問題.為得到方程(1)的穩(wěn)定數(shù)值解，本文采用了多重離散化光滑條件，所采用的方法更直觀、靈活.具體的理論方法見本文的第2部分和第3部分，第4部分的數(shù)值模擬驗(yàn)證本文所采用方法的可行性和有效性.

2 積分方程的離散化

下面將利用三次樣條插值函數(shù)對(duì)積分方程(1)進(jìn)行離散.首先，將積分區(qū)間[a，b]進(jìn)行等距剖分成n段，步長為h，每一個(gè)分點(diǎn)xi=a+ih，i=0，1，…，n，相應(yīng)地yj=a+jh，j=0，1，…，n，則在整個(gè)區(qū)間[a，b]上，積分方程(1)可以寫成如下形式

下面對(duì)未知函數(shù)f(x)應(yīng)用三次樣條插值函數(shù)S(x)進(jìn)行逼近，由于S(x)∈C2[a，b]，于是可設(shè)S(x)在xi處的一階導(dǎo)數(shù)值為mi，即S′(xi)=mi.在小區(qū)間[xi，xi+1]上，S(x)是一個(gè)不高于三次的多項(xiàng)式，記S(x)在區(qū)間[xi，xi+1]上的表達(dá)式為S3(x)，于是有

S3(xi)=f(xi)，S3(xi+1)=f(xi+1);S′3(xi)=mi，S′3(xi+1)=mi+1.

實(shí)際上，S3(x)為滿足以上條件的一個(gè)三次Hermit插值，則根據(jù)Hermit插值的形式得

S3(x)=hi(x)fi+hi+1fi+1+Hi(x)mi+Hi+1(x)mi+1，

(2)

其中，fi和fi+1分別為f(x)在xi和xi+1上函數(shù)值，hi(x)，hi+1(x)，Hi(x)，Hi+1(x)分別為三次Hermit插值的基函數(shù)，形式為

hi(x)=(h+2(x-xi))(x-xi+1)2/h3，hi+1(x)=(h+2(x-xi+1))(x-xi)2/h3，

將S3(x)的表達(dá)式(2)帶入方程(1)中可得

(3)

其中

將方程(3)寫成如下的矩陣向量的形式：

H(f0，f1，…，fn)T+H*(m0，m1，…，mn)T=(g0，g1，…，gn)T.

(4)

由S″(x)的表達(dá)式，可得

2m0+m1=3f(x0，x1)-h1M0/2，mn-1+2mn=3f(xn-1，xn)+hnMn/2，

其中f(x0，x1)表示x0，x1兩點(diǎn)處的一階差商，f(xn-1，xn)表示xn-1，xn兩點(diǎn)處的一階差商.當(dāng)M0，Mn值都為0時(shí)，稱為自然邊界條件，適合自然邊界條件的樣條稱為自然樣條.此時(shí)，方程(3)可以寫為

(5)

其中

將λi，μi帶入式(5)，并將差分格式f(xi-1，xi)展開，有

B(m0，m1，…，mn)T=K(f0，f1，…，fn)T，

(6)

其中

則

(m0，m1，…，mn)T=B-1K(f0，f1，…，fn)T，

將上式代入到式(4)，有

(H+H*B-1K)(f0，f1，…，fn)T=(g0，g1，…，gn)T.

(7)

不妨將(7)式簡記為

Af=g，

(8)

其中A=(H+H*B-1K)，f=(f0，f1，…，fn)T，g=(g0，g1，…，gn)T.方程(8)就是具有二階連續(xù)導(dǎo)數(shù)的積分方程的數(shù)值求解模型，不難發(fā)現(xiàn)，這是一類第一類算子方程.求解(8)式就可以得到積分方程(1)的近似解.

3 正則化策略

由于第一類Fredholm積分方程的不適定性，導(dǎo)致其離散化后得到的線性方程組(8)的解在Hadarm意義下是不適定的，因此為了得到方程(8)的適定解，一般采用Tikhonov正則化策略，即將對(duì)(8)的求解轉(zhuǎn)化成如下的無約束極小化問題

(9)

(10)

也可以采用單一罰項(xiàng)的約束.本文將在下一節(jié)對(duì)不同約束進(jìn)行對(duì)比.對(duì)公式(10)右端的三項(xiàng)采用離散近似，即

并將近似表達(dá)式帶入到(9)中，則重構(gòu)出來的光滑解實(shí)際上就是使下式達(dá)到最小值的f

(11)

其中λ1，λ2，λ3是三個(gè)正則化參數(shù).對(duì)(11)中的目標(biāo)泛函關(guān)于f求極小，并用矩陣形式表示，則有

AT(Af-g)+λ0H0f+λ1H1f+λ2H2f=0，

(12)

其中H0就是單位陣，Hi，i=1，2為加入的1階和2階導(dǎo)數(shù)約束所對(duì)應(yīng)的光滑矩陣，其形式分別如下：

求解方程(12)，整理可得

f=(ATA+λ0H0+λ1H1+λ2H2)-1ATg，

即為積分方程(1)的數(shù)值近似解.

4 數(shù)值模擬計(jì)算

通過數(shù)值模擬來檢驗(yàn)本文提出的方法的可行性與穩(wěn)定性.首先生成數(shù)據(jù)，即先給定積分方程的精確的解析形式f(x)，并計(jì)算出積分方程的右端函數(shù)g(y)，及其在不同的離散點(diǎn)上的值g(yi)，i=1，…，n.值得指出的是，本文中主要考慮精確解是光滑連續(xù)函數(shù)的情形.適當(dāng)?shù)貙?duì)數(shù)據(jù)加入人為的干擾誤差，以模擬實(shí)際數(shù)據(jù)中的測量誤差.在下面的例子中，對(duì)g(yi)加入了0～1%的高斯白噪聲.

以相對(duì)誤差RE=(‖f*-f‖)/‖f‖，和殘差RESIE=‖Af*-g‖為評(píng)價(jià)指標(biāo)，對(duì)比了多重約束與單一約束對(duì)重構(gòu)解的不同影響.同時(shí)，由于算法中有三個(gè)不同的正則化參數(shù)，為了找出最優(yōu)的參數(shù)，首先經(jīng)驗(yàn)地給出H2前的正則化參數(shù)λ2，在此基礎(chǔ)上，對(duì)H0和H1前面的參數(shù)設(shè)計(jì)了循環(huán)遍歷算法，以找出能使相對(duì)誤差達(dá)到極小的參數(shù)組合.

例1真解f(x)=sin(x)，核函數(shù)A(x，y)=exy，x∈[1，2]，不同約束條件下的重構(gòu)結(jié)果見圖1所示，其中單獨(dú)約束的情況正則化參數(shù)均為0.005，剖分的節(jié)點(diǎn)總數(shù)n=30. 混合約束中，二階導(dǎo)數(shù)H2前面的正則化參數(shù)為0.005，在此基礎(chǔ)上，對(duì)0階導(dǎo)數(shù)和1階導(dǎo)數(shù)前面的正則參數(shù)各循環(huán)遍歷10次，參數(shù)遍歷的范圍為初值為0.5，公比為0.1的等比數(shù)列.根據(jù)這100種(λ0，λ1，λ2)的組合，可以得到相應(yīng)的重構(gòu)解，并可以計(jì)算每一種重構(gòu)解的相對(duì)誤差.選擇相對(duì)誤差最小的，對(duì)應(yīng)的組合為

λ0=5×10-6，λ1=5×10-8，λ2=0.005.

根據(jù)圖1，不難看出混合約束的重構(gòu)效果是最好的，精確解與重構(gòu)解幾乎是重合的.在采用單獨(dú)約束的方法中，二階導(dǎo)數(shù)由于具有較強(qiáng)的光滑性，因此其重構(gòu)效果相較圖1(c)和圖1(d)要好.為了更好地量化對(duì)比的結(jié)果，列出了不同方法下重構(gòu)解的相對(duì)誤差和殘差，見表1.通過表1，可以看出多重約束方法具有最低的相對(duì)誤差，這也說明了采用多重約束確實(shí)能提高重構(gòu)解的精度.

圖1 不同約束下精確解與重構(gòu)解的對(duì)比

表1 不同約束下重構(gòu)解的相對(duì)誤差和殘差

例2真解f(x)=sin(5x)，核函數(shù)A(x，y)=exy，x∈[1，2]，通過例1可以發(fā)現(xiàn)混合約束與二階導(dǎo)數(shù)約束的重構(gòu)均優(yōu)于一階導(dǎo)數(shù)約束和0階導(dǎo)數(shù)約束.本例中僅對(duì)比混合約束與二階導(dǎo)數(shù)約束，結(jié)果見圖2(a)，其中二階導(dǎo)數(shù)約束前面的正則化參數(shù)為0.005，混合約束的正則化參數(shù)組合為λ0=5×10-6，λ1=5×10-6，λ2=0.005.兩種罰項(xiàng)約束下得到的相對(duì)誤差和殘差見表2的第二行.

表2 例2和例3不同精確解下重構(gòu)解的相對(duì)誤差和殘差

圖2 例2和例3的重構(gòu)結(jié)果

例3真解f(x)=ex，核函數(shù)同上.仍然對(duì)比混合約束與二階導(dǎo)數(shù)約束的重構(gòu)效果.結(jié)果見圖2-(b)，其中二階導(dǎo)數(shù)約束前面的正則化參數(shù)為0.005，混合約束的正則化參數(shù)組合為λ0=5×10-7，λ1=5×10-8，λ2=0.005.兩種罰項(xiàng)約束下得到的相對(duì)誤差和殘差見表2的第三行.

通過例2和例3，可以發(fā)現(xiàn)，當(dāng)精確解在重構(gòu)區(qū)間上，滿足單調(diào)性時(shí)，如例3，多重光滑約束和二階導(dǎo)數(shù)約束的結(jié)果幾乎無異，甚至二階導(dǎo)數(shù)的誤差量化結(jié)果要優(yōu)于多重光滑約束；但當(dāng)精確解在重構(gòu)區(qū)間上，雖然連續(xù)光滑，但不滿足單調(diào)性時(shí)，如例2，則多重光滑約束的重構(gòu)效果要優(yōu)于二階導(dǎo)數(shù)約束，特別是在相對(duì)誤差的量化上，多重光滑約束的相對(duì)誤差幾乎是二階導(dǎo)數(shù)約束相對(duì)誤差的2倍.

5 結(jié) 論

本文利用三次樣條插值函數(shù)對(duì)積分方程進(jìn)行離散，并對(duì)離散得到的線性方程組采用光滑約束，以得到穩(wěn)定解.得到結(jié)論如下：

(i) 當(dāng)精確解的光滑性較高時(shí)，采用三次樣條函數(shù)對(duì)其進(jìn)行近似符合對(duì)解的光滑性的要求，且方法簡單易行，轉(zhuǎn)化成無約束問題后直接求極小解線性方程組即可，不用進(jìn)行迭代求解；

(ii) 當(dāng)精確解的光滑性較高，且滿足一定的單調(diào)性時(shí)，采用多重光滑約束可以得到精度較高的重構(gòu)解；

(iii) 方法中的正則化參數(shù)的選取，要視對(duì)解的要求而定，針對(duì)不同要求，選取不同的光滑參數(shù)，如果解的光滑性要求較高，則應(yīng)以高階導(dǎo)數(shù)約束為主，適當(dāng)減小其他低階導(dǎo)數(shù)約束前面的正則化參數(shù).

致謝作者非常感謝相關(guān)文獻(xiàn)對(duì)本文的啟發(fā)以及審稿專家提出的寶貴意見.