Fredholm積分方程的數(shù)值解

1、第18卷第2期核聚變與等離子體物理V ol.18,N o.2 1998年6月N uclear Fusion and Pla sma Physics J une1998Fredholm積分方程的數(shù)值解徐文斌董家齊(核工業(yè)西南物理研究院,成都610041積分方程是泛函分析的一個(gè)重要分支,它是研究數(shù)學(xué)及其它學(xué)科,如偏微分方程邊值問題和各種物理問題的重要數(shù)學(xué)方法。為此,我們研制了求解Fredho lm積分方程的代碼IT GM O,并給出了該代碼實(shí)際應(yīng)用的例子。關(guān)鍵詞Fredho lm積分方程Z i模數(shù)值解1引言在等離子體物理的計(jì)算中,經(jīng)常需求解Fredholm積分方程,如求解Z i模的不穩(wěn)定性1,4和

2、求解微撕裂模的不穩(wěn)定性2時(shí)就會(huì)碰到這樣的問題。為此我們編制了求解Fredholm積分方程的代碼,以便研究和分析各種有關(guān)的物理問題。積分方程的一般形式為:T(sx(s=y(s+b a k(s,tx(td t(1式中T(s、y(s和k(s,t是已知函數(shù),、a、b是變數(shù)。s和t可取a,b區(qū)間上的一切值。稱為積分方程的參數(shù),x(s稱為積分方程的解函數(shù),k(s,t稱為積分方程的核,y(s項(xiàng)稱為自由項(xiàng)。如果方程(1的未知函數(shù)是一次的,則稱方程(1為線性積分方程。如果y(s0,就稱方程(1為齊次的積分方程,否則稱為非齊次積分方程。當(dāng)y(s0時(shí),又可稱為積分方程(1的本征值,x(s為本征函數(shù)。一般將Fredh

3、olm方程分為三類:第一類方程為:b a k(s,tx(td t=y(s(2第二類方程為:x(s=y(s+b a k(s,tx(td t(3第三類方程為:T(sx(s=y(s+b a k(s,tx(td t(4當(dāng)?shù)谌惙匠痰腡(s在a,b上是正函數(shù)時(shí),此類方程便可以化成:國(guó)家自然科學(xué)基金(19575014資助26核聚變與等離子體物理第18卷T(sx(s=y(sT(s+b a k(s,tT(sT(tT(tx(td t(5因此它變成含有x(s=T(sx(s的第二類Fredholm方程。另外,在線性方程(3中,如果y(s0,則方程(3的通解可寫成x=x0+Z,其中x0是相應(yīng)齊次方程的解,Z為方程(3

4、的一特解。這就說(shuō)明了只要求出齊次方程的解,就很容易得到非齊次方程(3的解。因此為了便于討論,我們?cè)谶@里取方程(3中的y(s0。下面我們重點(diǎn)來(lái)研究線性齊次第二類Fredholm 積分方程的數(shù)值解。2第二類Fredhol m積分方程的解法當(dāng)y(s0時(shí),令V=1,則方程(3可寫為:k(s,tx(td t=V x(s(6將方程(6寫成算子的形式,則有:Kx=V x(7式中K是積分算子,V是本征值,x為本征函數(shù)。在L2空間中按Releigh-Ritz方法定義一內(nèi)積函數(shù):F(u=(Ku,u(u,u(8顯然,當(dāng)x是方程(7的解函數(shù)時(shí),有:F(x=(Kx,x (x,x=V這說(shuō)明這種Rele ig h函數(shù)等價(jià)于

5、積分算子K。當(dāng)函數(shù)u與真實(shí)的本征函數(shù)x有一小的誤差e時(shí),即:u=x+e(9則有:F(u=F(x+e=(Kx,u+(Ke,u(u,u=(V x,u+(Ke,x+(Ke,e(u,u=V(u,u-(x,e-(e,e+V(x,e+(Ke,e(u,u=V+e,(K-V Ie(u,u(10方程(10給出了關(guān)于V的Releig h商。由方程(10可知,當(dāng)函數(shù)x有一小的誤差時(shí),本征值的精度為e2。不僅如此,如果V max是函數(shù)F的最大本征值,那么對(duì)任意的e,有e,(K-V Ie0。因此,對(duì)任意的向量x,有:F(xV max(11第2期徐文斌等:Fredholm 積分方程的數(shù)值解27為了求解積分方程(6,可選擇

6、n 個(gè)函數(shù)h 1,h 2,h n ,用它們的線性組合來(lái)逼近解向量x ,即:x (s x n (s =ni =1a i h i(s (12其中n 個(gè)系數(shù)a =(a 1,a 2,a n T 可用最小二乘法來(lái)確定。由式(8可知:F (x n =(Kx n ,x n (x n ,x n =ni =1nj =1a i a j h i(s K (s ,t h j(t d s d tni =1nj =1a i a jh i (s h j (t d s d t(13為了方便起見,可把式(13寫成下式:F (x n =a TLaa TMa(14式中L ij =h i (s K (s ,t h j (t d s

7、d tM ij=h i(s h j(t d s d t (i ,j =1,2,n (15通過對(duì)F (x n 取極值,就可得到使e n =x -x n 為最小的a 向量。令 F / a i =0,可得到:(a TMa La -(a TLa Ma (a TM a =0從而有:La -V n Ma =0(16當(dāng)h i (s (i =1,2,n 取成歸一化的正交函數(shù)時(shí),M ij 就變成了單位矩陣I ,這時(shí)式(16變成:La =Vn a (17式中V n =F (x n 是相應(yīng)于積分算子K 的近似矩陣L 的本征值。由式(11可知,式(17的最大本征值V maxn =max V 1n ,V 2n ,V n

8、n 恒小于或等于式(7中的真實(shí)最大本征值V max,這就是說(shuō)式(17的最小本征值minn =1/V maxn 恒大于或等于式(3中的最小本征值min。因此,只要求出式(17的最大本征值V max n 和相應(yīng)的本征向量x n ,就找到了原方程(6的解。為了得到式(17的最大本征值V maxn 和本征向量x n ,可采用改進(jìn)的冪法3進(jìn)行求解。冪法是從任取的初始向量x(0出發(fā),用矩陣逐次乘這個(gè)向量,即:x (0,x(1=Lx(0,x(2=Lx(1,x(i =Lx(i -1(18由式(12可知,x 可表成n 個(gè)特性向量的線性組合,因此有:x (i =(i 1a 1h 1+a 2(21i h 2+a 3

9、(31i h 3+a n (n 1i h n (19這說(shuō)明初始向量x (0在特征向量x j 上的分量按(j1k 的速度收斂到零。假設(shè)1和2分別為最大和次大的本征值,那么當(dāng)2/1接近于1時(shí),收斂速度就會(huì)很慢。為了加快收斂速度,用相 28核聚變與等離子體物理第18卷鄰三次近似向量的組合來(lái)實(shí)現(xiàn)W 2加速。其辦法如下:假設(shè)迭代過程已進(jìn)行到h 3-h n 上的分量可以忽略的程度,那么y (i =x (i /max (x (i 將與h 1+W h 2相差一個(gè)常數(shù)因子。由于W 很小,于是有max (h 1+W h 2=1+W p ,這里p 為h 2的分量中與h 1按模最大的相應(yīng)者。此外,以后的近似向y(i

10、+1中按模最大的分量的位置亦不再改變,這樣就可以把y (i 、y(i +1、y(i +2寫為:(h 1+W h 2/(1+W p ,(1h 1+W 2h 2/(1+W p 2,(21h 1+W 22h 2/(21+W p 22現(xiàn)考慮如下向量Z =(Z 1,Z 2,Z n T ,其中Z j =y (i j y (i +2j -(y (i +1j 2/(y (i j -2y (i +1j +y (i +2j =y (i j -y (i j -y (i +1j 2/y (i j -2y (i +1j +y (1+2j (j =1,2,n (20其中y (i j 為向量y (i 的第j 個(gè)分量。把前面

11、關(guān)于y (i 、y (i +1、y(i +2的表達(dá)式代入上式,并以h (i j 表示h i 的第j 個(gè)分量,則可以得到:Z j =h (1j-W 2(21ph (2j /1-W 2p 2(212(21這就是說(shuō),Z =h 1+O (W 2,因而對(duì)h 1來(lái)說(shuō),Z 是較y (k +2更為精確的近似。采用這種方法來(lái)求式(7的本征值和本征函數(shù),實(shí)際計(jì)算中一般只要迭代10步就可找到收斂的解。3計(jì)算例子與具體解法3.1作為例子的方程為了研究電子溫度梯度驅(qū)動(dòng)的短波模,Dong 1導(dǎo)出了在k 空間求解擾動(dòng)靜電勢(shì)的積分方程:T (k h (k +K (k ,k h (k d k =0(22這顯然是齊次的第三類Fr

12、edholm 積分方程,其中T (k =f (1+1f z 1+0(k i +k 2i Zi i 2f z0(k i -1(k i y (k =0K (k ,k =12ex pi x (k -k (2x L m i /m e 22(e - e +m e m ie - i d x(23這里相應(yīng)于電子和離子的電導(dǎo)為:e - e =2e 1-(1/z 1+e Z (e -Z e /z 2e +(2e -12e Z (e (24e - i =2i (1+1f z1+i Z (i 0(k i ,k i +Z i f z 2i +(2i -12i Z (i 0(k i ,k i +Z i f z 1+i

13、Z (i k i k i 2f 1(k i ,k i -14f(k 2i +k 2i 0(k i ,k i (25第2期徐文斌等:Fredholm 積分方程的數(shù)值解29其中,Z (e 和Z (i 分別為電子和離子的色散函數(shù),f =T e T i ,z =k k *ne ,k *ne =K y T e C eBL n,L n =(-1n d n d x -1,j (k i =I j (k i 2f e x p (-k 2i2f,k =k r +i V ,j (k i ,k i=I j (k i ,kie x p (-k 2i +k 2i 2f,k 2i =k 2y +k 2,k 2i =k 2y

14、 +k 2,I j (k i 和I j (k i ,k i 是第j 階的修正的貝塞爾函數(shù),e =z 2x Lm i /m e,i =f z 2x L ,L =L n L s,Z i =d(ln T i d(ln n ,Z e =d(ln T e d(ln n 。在這些方程中所有的長(zhǎng)度都用d e =2T e /m i K 2i 1/2進(jìn)行了歸一化,且假定k d e 1。方程(17支配了短波模h (k 的性質(zhì)。我們的任務(wù)就是求解h (k 。這里x 變量的關(guān)系與k 變量的關(guān)系為傅里葉的變換關(guān)系,即有:p (x =12p (k ex p(i x R d k(263.2具體解法從上面對(duì)方程的描述可知,該

15、方程的T (k 和K (k ,k 項(xiàng)很復(fù)雜,如采用變換的方法將式(22變?yōu)榈诙怓redholm 積分方程來(lái)求解,則K (k ,k 項(xiàng)會(huì)變得更加復(fù)雜。為此我們不對(duì)式(22進(jìn)行變換,而采用直接積分的離散方法進(jìn)行求解。令a 0,a 1,a n 為k 空間上的等分點(diǎn),h 為任一小區(qū)間a i ,a i +1的長(zhǎng)度(i =0,1,n -1,則在小區(qū)間a i ,a i +1上,式(22可寫為:ai +1aiT (k h (k d k +a i +1aiK (k ,k h(k d k d k =0(27采用矩形積分法:ai +1a iT (k h (k d k =T (a i h (a i ha i +1a

16、iK (k ,k h (k d k d k =h K (a i ,k h (k d k =h 2n -1j =0K (a i ,a j h (a j因此在區(qū)間a i ,a i +1上,原方程化為:h T (a i h (a i +h 2K i H =0(28其中,K i =K (a i ,a 0,K (a i ,a n -1,H =h (a 0,h (a 1,h (a n -1T 。對(duì)所有的i =0,1,n -1,方程可整理為:hT (a 00T (a 10T (a n -1+h 2K 0K 1K n -1H =A ij H =0(29當(dāng)A i j 是一個(gè)大的矩陣時(shí),為了方便地求解式(29,我

17、們可替代地用W 2加速方法解下式:30 核聚變與等離子體物理 第 18 卷 det ( Aij - _W ij = 0 _ 就是矩陣 Ai j 的本征值 , 當(dāng) _ 趨于零時(shí) , 就求得了方程 ( 29 的解 z = k / k 。 * ne ( 30 3. 3 計(jì)算結(jié)果與討論 V k r 圖1 示出了用我們編制的 ITGM O代碼計(jì)算出的 Z i 模的增長(zhǎng)率 * 和實(shí)頻率 * 隨 Z i 值的變 k ne k ne 化。 圖1 Z i 模的增長(zhǎng)率和實(shí)頻率隨 Z i 的變化 - 1 = 2. 5× 10- 2。 k y = 1, Z e = 2;f= 1; L n L s 圖 2 示

18、出了典型的 Z i 模在 k 空間和 x 空間的本征函數(shù)。 圖 2 Z i 模在 k 空間和 x 空間的典型曲線 - 1 k y = 0. 35; f = 1; Z 。 i = 4; L n L s = 0. 15 總之 , 我們研制的求解積分方程的代碼是很有效的計(jì)算工具。在等離子體物理中 ,該代碼 擴(kuò)展后還可用于求解撕裂模、環(huán)形條件下的 Z i 模和有雜質(zhì)情況的 Z i 模等許多實(shí)際的物理問題。 (下轉(zhuǎn)第 51 頁(yè) 第 2期 祁蘭英等 : 激光超熱電子與受激 Raman 散射特性的實(shí)驗(yàn)觀測(cè) 51 T he energy and spect rum characteristic of hot

19、electrons and stimulated Raman backward sca ttering are experim entally studied. A eff ective w ay to suppress ho t elect rons has been f ound. Key words Gold-disk targ et Hohlraum target Hydrocarbon f oil t arget Hot electron Stimulated Ram an scat tering (上接第 30 頁(yè) 參考文獻(xiàn) 1 Dong J Q, Gu zdar P N , Le

20、e Y C. Fini te Bet a Ef fects on Ion Tem perat ure G radi ent Driv en M odes. Phys. Fluids, 1987, 30( 9: 2694. 2 Farengo R, Lee Y C, G uzdar P N. An Elect romagn etic Int egral Equation : Application to M icrot earing M odes. Ph ys. Fluids , 1983, 26( 12: 3515. 3 馮康 . 數(shù)值計(jì)算方法 . 北京: 國(guó)防工業(yè)出版社 , 1978. 402. 4 Dong J