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積分方程是泛函分析的一個重要分支,它是研究數(shù)學及其它學科,如偏微分方程邊值問題和各種物理問題的重要數(shù)學方法。為此,我們研制了求解Fredho lm積分方程的代碼IT GM O,并給出了該代碼實際應用的例子。關鍵詞Fredho lm積分方程Z i模數(shù)值解1引言在等離子體物理的計算中,經常需求解Fredholm積分方程,如求解Z i模的不穩(wěn)定性1,4和

2、求解微撕裂模的不穩(wěn)定性2時就會碰到這樣的問題。為此我們編制了求解Fredholm積分方程的代碼,以便研究和分析各種有關的物理問題。積分方程的一般形式為: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項稱為自由項。如果方程(1的未知函數(shù)是一次的,則稱方程(1為線性積分方程。如果y(s0,就稱方程(1為齊次的積分方程,否則稱為非齊次積分方程。當y(s0時,又可稱為積分方程(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當?shù)谌惙匠痰腡(s在a,b上是正函數(shù)時,此類方程便可以化成:國家自然科學基金(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是相應齊次方程的解,Z為方程(3

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

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

6、n 個函數(shù)h 1,h 2,h n ,用它們的線性組合來逼近解向量x ,即:x (s x n (s =ni =1a i h i(s (12其中n 個系數(shù)a =(a 1,a 2,a n T 可用最小二乘法來確定。由式(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通過對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當h i (s (i =1,2,n 取成歸一化的正交函數(shù)時,M ij 就變成了單位矩陣I ,這時式(16變成:La =Vn a (17式中V n =F (x n 是相應于積分算子K 的近似矩陣L 的本征值。由式(11可知,式(17的最大本征值V maxn =max V 1n ,V 2n ,V n

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

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

11、關于y (i 、y (i +1、y(i +2的表達式代入上式,并以h (i j 表示h i 的第j 個分量,則可以得到:Z j =h (1j-W 2(21ph (2j /1-W 2p 2(212(21這就是說,Z =h 1+O (W 2,因而對h 1來說,Z 是較y (k +2更為精確的近似。采用這種方法來求式(7的本征值和本征函數(shù),實際計算中一般只要迭代10步就可找到收斂的解。3計算例子與具體解法3.1作為例子的方程為了研究電子溫度梯度驅動的短波模,Dong 1導出了在k 空間求解擾動靜電勢的積分方程: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這里相應于電子和離子的電導為: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 。在這些方程中所有的長度都用d e =2T e /m i K 2i 1/2進行了歸一化,且假定k d e 1。方程(17支配了短波模h (k 的性質。我們的任務就是求解h (k 。這里x 變量的關系與k 變量的關系為傅里葉的變換關系,即有:p (x =12p (k ex p(i x R d k(263.2具體解法從上面對方程的描述可知,該

15、方程的T (k 和K (k ,k 項很復雜,如采用變換的方法將式(22變?yōu)榈诙怓redholm 積分方程來求解,則K (k ,k 項會變得更加復雜。為此我們不對式(22進行變換,而采用直接積分的離散方法進行求解。令a 0,a 1,a n 為k 空間上的等分點,h 為任一小區(qū)間a i ,a i +1的長度(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 。對所有的i =0,1,n -1,方程可整理為:hT (a 00T (a 10T (a n -1+h 2K 0K 1K n -1H =A ij H =0(29當A i j 是一個大的矩陣時,為了方便地求解式(29,我

17、們可替代地用W 2加速方法解下式:30 核聚變與等離子體物理 第 18 卷 det ( Aij - _W ij = 0 _ 就是矩陣 Ai j 的本征值 , 當 _ 趨于零時 , 就求得了方程 ( 29 的解 z = k / k 。 * ne ( 30 3. 3 計算結果與討論 V k r 圖1 示出了用我們編制的 ITGM O代碼計算出的 Z i 模的增長率 * 和實頻率 * 隨 Z i 值的變 k ne k ne 化。 圖1 Z i 模的增長率和實頻率隨 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 總之 , 我們研制的求解積分方程的代碼是很有效的計算工具。在等離子體物理中 ,該代碼 擴展后還可用于求解撕裂模、環(huán)形條件下的 Z i 模和有雜質情況的 Z i 模等許多實際的物理問題。 (下轉第 51 頁 第 2期 祁蘭英等 : 激光超熱電子與受激 Raman 散射特性的實驗觀測 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 頁 參考文獻 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ù)值計算方法 . 北京: 國防工業(yè)出版社 , 1978. 402. 4 Dong J