Fredholm积分方程的数值解

1、第18卷第2期核聚变与等离子体物理V ol.18,N o.2 1998年6月N uclear Fusion and Pla sma Physics J une1998Fredholm积分方程的数值解徐文斌董家齐(核工业西南物理研究院,成都610041积分方程是泛函分析的一个重要分支,它是研究数学及其它学科,如偏微分方程边值问题和各种物理问题的重要数学方法。为此,我们研制了求解Fredho lm积分方程的代码IT GM O,并给出了该代码实际应用的例子。关键词Fredho lm积分方程Z i模数值解1引言在等离子体物理的计算中,经常需求解Fredholm积分方程,如求解Z i模的不稳定性1,4和

2、求解微撕裂模的不稳定性2时就会碰到这样的问题。为此我们编制了求解Fredholm积分方程的代码,以便研究和分析各种有关的物理问题。积分方程的一般形式为:T(sx(s=y(s+b a k(s,tx(td t(1式中T(s、y(s和k(s,t是已知函数,、a、b是变数。s和t可取a,b区间上的一切值。称为积分方程的参数,x(s称为积分方程的解函数,k(s,t称为积分方程的核,y(s项称为自由项。如果方程(1的未知函数是一次的,则称方程(1为线性积分方程。如果y(s0,就称方程(1为齐次的积分方程,否则称为非齐次积分方程。当y(s0时,又可称为积分方程(1的本征值,x(s为本征函数。一般将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当第三类方程的T(s在a,b上是正函数时,此类方程便可以化成:国家自然科学基金(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 积分方程的数值解。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为本征函数。在L2空间中按Releigh-Ritz方法定义一内积函数:F(u=(Ku,u(u,u(8显然,当x是方程(7的解函数时,有:F(x=(Kx,x (x,x=V这说明这种Rele ig h函数等价于

5、积分算子K。当函数u与真实的本征函数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可知,当函数x有一小的误差时,本征值的精度为e2。不仅如此,如果V max是函数F的最大本征值,那么对任意的e,有e,(K-V Ie0。因此,对任意的向量x,有:F(xV max(11第2期徐文斌等:Fredholm 积分方程的数值解27为了求解积分方程(6,可选择

6、n 个函数h 1,h 2,h n ,用它们的线性组合来逼近解向量x ,即:x (s x n (s =ni =1a i h i(s (12其中n 个系数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 取成归一化的正交函数时,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出发,用矩阵逐次乘这个向量,即: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卷邻三次近似向量的组合来实现W 2加速。其办法如下:假设迭代过程已进行到h 3-h n 上的分量可以忽略的程度,那么y (i =x (i /max (x (i 将与h 1+W h 2相差一个常数因子。由于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现考虑如下向量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的本征值和本征函数,实际计算中一般只要迭代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 积分方程的数值解29其中,Z (e 和Z (i 分别为电子和离子的色散函数,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 阶的修正的贝塞尔函数,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变为第二类Fredholm 积分方程来求解,则K (k ,k 项会变得更加复杂。为此我们不对式(22进行变换,而采用直接积分的离散方法进行求解。令a 0,a 1,a n 为k 空间上的等分点,h 为任一小区间a i ,a i +1的长度(i =0,1,n -1,则在小区间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因此在区间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 空间的本征函数。 图 2 Z i 模在 k 空间和 x 空间的典型曲线 - 1 k y = 0. 35; f = 1; Z 。 i = 4; L n L s = 0. 15 总之 , 我们研制的求解积分方程的代码是很有效的计算工具。在等离子体物理中 ,该代码 扩展后还可用于求解撕裂模、环形条件下的 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 冯康 . 数值计算方法 . 北京: 国防工业出版社 , 1978. 402. 4 Dong J