fredholm,离散积分方程

1、1 第一类 Fredholm 积分方程，具有形式如下：b(1)k( x, s) y(s)ds f (x) , a x ba其中核函数 K ( x, s) 和自由项f ( x) 为已知函数，y(s) 是未知函数。此类积分方程虽然形式简单， 但其求解却比较困难， 所以这类方程在下文将做详细介绍。2 第二类 Fredholm 积分方程，具有如下的形式：y( x)b(2)k (x, s) y( s) ds f ( x) , a x ba离散积分方程的数值方法有很多种，比如可以用复化梯形公式、复化辛普森公式等，这里我们利用复化梯形公式来进行离散。一、复化梯形公式离散过程如下：bh f ( a)nf (x

2、)dx2f ( xk ) f (b)a2k1下面具体给出复化梯形公式对第二类积分方程的一般离散过程。bs1s2k( x, s) y(s)dsk (x, s) y(s)dsk( x, s) y(s)dsas0s1sisnk( x, s) y( s)dsk( x, s) y(s)dssi1sn 1h k (x, s0 ) y(s0 )k (x, s1 ) y(s1 )2hk( x, sn ) y( sn )k (x, sn 1 ) y( sn 1 )2h 1 k( x, s0 ) y(s0 )k( x, s1 ) y( s1 )2 h2 (b a) k(x,) f ()12h k (x, si 1

3、 ) y(si1 )k (x, si ) y(si )2h 2 (b a) f ()k (x,12k( x, si ) y( si )1 k (x, sn ) y(sn )2g (x) a,b最后对变量 x 进行离散，将区间 a, b 等分为 n 份，步长为 hb a ，n同时忽略积分公式误差项：g( xi )y(xi )h 1 k (xi , s0 ) y(s0 )k ( xi , s1 ) y( s1 )k(xi , si ) y(si )21 k ( xi , sn ) y(sn )2其中 i0,1,2,n得到线性方程组Af ng n其中 fn y( s0 ), y(s1 ), y(s2

4、 ),y(sn ) ， gn g( x0 ), g( x1 ), g (x2 ), g (xn )1hk( x0 , y1 )1 hk (x0 , y0 )212hk (x0 , yn )A1hk( x1 , y0 )1hk (x1 , y1 )212hk( x1 , yn )1hk( xn , y1 )h k( xn , y0 )211 hk( xn , yn )2再对上述方程进行数值求解，即可。1(1例：求解积分方程 y( x)x) y( s)ds01s4x 35x 22x58( x1) 2，其解析解为y( x)(1x) 2代码如下：functionK=K(x,y)K = 1/(1+y)

5、- x;functionw1=fun1(x)w1=1./(1+x).*(1+x);functionf=f(x)functionw5=fww(a,b,n)%第一类 fredholm 方程解的程序%w5=w1,w2,w3,w4,各列分别表示真解、数值解、最小二乘解、正则解%a,b表示积分区间a,b%n表示将区间 n等分%m表示正则参数的取值h=(b-a)/n;x=a:h:b;y=a:h:b;A=zeros(n+1,n+1);%初始化矩阵 A为 n+1阶零矩阵g=zeros(n+1,1);%初始化列向量g为 n+1维零向量w1=zeros(n+1,1);%初始化列向量w1为 n+1维零向量fori=

6、1:n+1for j=1:n+1A(i,j)=K(x(i),y(j);endg(i)=f(x(i);w1(i)=(fun1(x(i);%计算方程的真解endA(:,1)=A(:,1)/2;A(:,n+1)=A(:,n+1)/2;A=h*A;A=eye(n+1,n+1)-A;w2=Ag; %得到的数值解aa=norm(w1-w2)/norm(w1);%相对误差bb=norm(w1-w2);%绝对误差cc=w1 w2;plot(x,w1,'b+') %真解holdonplot(x,w2,'r*') %数值解%axis(0 1 -100 100);%设置坐标轴titl

7、e(' 数值解与真解的比较' ); %加图形标题xlabel(' 变量 y' ); %加 x轴说明ylabel('y 对应的解 ' ); %加 y轴说明运行结果：>> fww(0,1,50) aa =%相对误差bb =0.000693865370887685 %绝对误差数值解与真解的比较1.21.110.90.8解的应 0.7对y0.60.50.40.30.200.10.20.30.40.50.60.70.80.91变 量 y二、辛普森公式离散过程如下：bha bf (x)dx f (a) 4 f () f (b)a62下面给出复化梯

8、形公式对第二类积分方程的一般离散过程。由于辛普森公式中取到中点的值，所以我们在区间a, b 上取 2n1个点。bs1s2skk( x, s) y(s)dsk (x, s) y(s) dsk ( x, s) y(s)dsk(x, s) y(s)dsas0s1sk 1h k( x, s0 ) y( s0 ) 4k ( x, s0s1 ) y( s0s1 ) k ( x, s1 ) y( s1 )622hs1s2s1s2) k ( x, s2 ) y( s2 ) k( x, s1 ) y( s2 ) 4k( x,2) y(26h k( x, sk 164(ba)h) y( sk 1 )f (4) (

9、 )4k (x, sk 1sk ) y( sk 1 sk ) k (x, sk ) y(sk )22h k( x, s0 ) y( s0 )4k (x, s1 ) y(s1 )2k( x, s2 ) y( s2 )4k( x, s2 n 1 ) y( s2n 1 ) k (x, s2n ) y( s2 n )(b a)h 4f ( 4) ( )2880g( x) a, b最后对变量x 进行离 散，将区 间 a, b 等分 为 2n 份，步 长为h b a ，同时忽略积分公式误差项：2ng(xi )y( xi )h k (xi , s0 ) y( s0 ) 4k (xi , s1 ) y(s1

10、) 2k( xi , s2 ) y( s2 )64k( xi , s2n 1 ) y(s2 n 1 )k ( xi , s2n ) y(s2 n )其中 i0,1,2,2n得到线性方程组 Af 2ng2n其中 f2 n y(s0 ), y( s1 ), y( s2 ), y( s2 n ) ， g2n g (x0 ), g( x1 ), g( x2 ), g (x2n )k( x0 , y0 )4k( x0 , y1 )2k( x0 , y2 )k( x0 , y2n )A Ih k( x1 , y0 )4k (x1, y1 )2k( x1 , y2 )k( x1 , y2n )6k( x2

11、n , y0 ) 4k( x2 n , y1 ) 2k ( x2 n , y2 )k (x2n , y2 n )11 hk( x0 , y0 )2 hk ( x0 , y1 )1 hk( x0 , y2 )1 hk (x0 , y2n )63361 hk ( x1 , y0 )12 hk (x1 , y1 )1 hk( x1 , y2 )1 hk ( x1 , y2n )63361 hk (x2n , y0 )2 hk ( x2 n , y1 )1 hk (x2n , y2 )11 hk (x2n , y2 n )6336再对上述方程进行数值求解，即可。例：求解积分方程 y( x)( 1x)

12、y(s)ds4x35x22x 5 ，其解101s8( x1) 2析解为 y(x)(1x) 2代码如下：functionv = knl(x,t)v = 1/(1+t) - x;functionf = fnc(x)functiony = inteqn(t, kernel, fun, coef)% Inputs% t evaluation points of the quadrature rule% kernel kernel function K% fun function f% coef quadrature rule coefficients% Output% y discrete solut

13、ion values at tn = length(t);f = feval(fun, t);%forj=1:nfori = 1:nK(j,i) = feval(kernel, t(j), t(i);endend% A = eye(n) - K*diag(coef); for j=1:nA(:,j) = -coef(j)*K(:,j);A(j,j) = 1.0 + A(j,j);endy = Af;k = input('Enter number of pannels: ');x = linspace(0,1,k+1); x = x'% evenly spaced kno

14、ts% Note: this x can be replaced by any partition of 0,1y = zeros(length(x),1);% discrete approximation at x% use Simpson sn = 2*k + 1;% number of points in Simpsons rulecoef = zeros(n,1);% coeficients in Simpsontrulet = zeros(n,1);% points in Simpsons rule% generate Simpsons rule coefficients and e

15、valuation pointsfori=1:kcoef(2*i-1:2*i+1) = coef(2*i-1:2*i+1) + 1; 4; 1;t(2*i-1) = x(i);t(2*i) = (x(i) + x(i+1)/2;endt(n) = x(k+1);coef = coef/(6*k);% solve the integral equationyt = inteqn(t, knl, fnc, coef);% discrete approximation to y(x) at the partition xy = yt(1:2:n);% check resultsyexact = 1.

16、/(1+x).*(1+x);aa=norm(y-yexact)/norm(y)%相对误差aa=norm(y - yexact)%绝对误差cc=y yexactplot(x,y,'b+') %真解holdonplot(x,yexact,'r*') %数值解%axis(0 1 -100 100);%设置坐标轴title(' 数值解与真解的比较' ); %加图形标题xlabel(' 变量 y' ); %加 x轴说明ylabel('y 对应的解 ' ); %加 y轴说明运行结果：Enter number of pannel

17、s: 50aa =6.9309212581323e-009aa =数值解与真解的比较1.21.110.9解0.8的0.7应对y0.60.50.40.30.200.10.20.30.40.50.60.70.80.91变 量 y三、高斯勒让德离散过程如下：关于定积分1 1,1, (x) 1，则关于权函数( x) f ( x)dx ，如果 a, b-1( x) 1 正交多项式就是1d m ( x21) nLn (x)dx n2 n n!这时 Gauss型积分公式的节点就取为上述多项式的零点，相应的Gauss型积分公式为1nf (x)dxAk f ( xk )-11k下面给出高斯公式对第二类积分方程的

18、一般离散过程。在 -1,1上 Fredholm 的积分公式为：1nAj k( x, sj ) y( s j)k( x, s) y(s)ds-1j 1g( x)第二类 Fredholm 积分方程可以化为：ng( xi )y( xi )Aj k( xi , s j ) y(sj )j11 A1k ( x1 , s1 )A2 k( x1 , s2 )An k( x1 , sn )y( s1 )g( x1 )即A1k(x2 , s1 )1 A2 k (x2 , s2 )An k( x2 , sn )y( s2 )g( x2 )A1k (xn , s1 )A2 k( xn , s2 )1 An k (x

19、n , sn )y(sn )g( xn )高斯勒让德型积分公式的积分区间为1,1 ，而对于一般的区间bf (x)dx 需要作变量替换 xba ba t 得到 a, b 上的积分a22bb a 1b abaf ( x)dx1f (t )dta222例：求积分方程111f (s)dsarctan(1x) arctan(1 s)y(x)11 ( x s)2其解析解 y( x) 1functionx,w = gauss(N)beta = .5./sqrt(1-(2*(1:N-1).(-2);T = diag(beta,1) + diag(beta,-1);V,D = eig(T);x = diag(D

20、) ;x,i = sort(x);w = 2*V(1,i).2;N=input('N=?' )ker= '(1/pi)*(1+(ss-tt).2).(-1)'s,w=gauss(N);t=s;ss,tt=meshgrid(s,t)ss=ss'tt=tt'K=eval(ker)W=diag(w);A=eye(N)+K*W;g=1+(1/pi)*(atan(1-s)+atan(1+s);ww=cond(K*W)u=Ag%数值解plot(s,u);运行结果：N=?5N =5u =0.9999835042915411.000043480651920.99

21、98914375874491.000043480651920.999983504291541四、克伦肖柯蒂斯（Clenshaw-curtis）离散过程如下：五、高斯罗巴托（ Gauss-Lobatto）离散：Gauss-Lobatto求积公式的表达式如下：12n 1f ( x)dx f (1) f ( 1)Aj f (x j ) Rn f -1n(n 1)j2Gauss-Lobatto求积公式的系数和余项分别为：Aj2n( n 1) Pn ( x j ) 2Rn f n 3 ( n 1)22 n 1( n 1)! 4f ( 2n ) ( ),( 1,1)(2n 1)( 2n)! 3其中， x

22、j 为 Pn (x) 的零点； Pn ( x) 为 n 次Legendre(勒让德 )多项式六、伽辽金（ Galerkin ）法离散：设i (x) 为 L2 a, b 空间内的一个完备正交系， k( x)L2 a,b ，则当 N充分大时，有nN ( x)kii( x)i0其中N ( x) 为( x) 的逼近函数。将上式带入 y( x)bk( x, s) y(s)ds g( x) 得annbki i ( x)kik (x, s) i (s)ds g (x)i 0i 0a两边分别对j (x) 求内积，得nbbki (ai ( x)k(x, s)i ( s),j (x) ( g( x), j (x)

23、i 0a即：nbb bbi0ai (x)j ( x)dxaak (x, s)i (s)j ( x)dsdxkiag (x)j (x)dxk0b0k1b1得：aijk2b2k nbn其中 aijbb bai ( x)j ( x)dxk (x, s) i ( s)j ( x) dsdxa abjbg ( x)j ( x)dxa例：求解第二类积分方程 y( x)xx y y(s)ds1 (e xe 3 x ) , 0 x 1 ，e02其解析解为 y(x) e x代码如下：functionw=obj(x,y)w=exp(-x-y);functionw=obj1(x)w=(exp(-x)+exp(-3*

24、x)/2;functionw1=phi_xk(x,k)ifk=0w1=ones(size(x);elseifk=1w1=x;elseifk=2w1=2*x.2-1;elseifk=3w1=4*x.3-3*x;elseifk=4w1=8*x.4-8*x.2+1;elseifk=5w1=16*x.5-20*x.3+5*x;elseifk=6w1=32*x.6-48*x.4+18*x.2-1;elseifk=7w1=64*x.7-112*x.5+56*x.3-7*x;elseifk=8w1=128*x.8-256*x.6+160*x.4-32*x.2+1;elseifk=9w1=256*x.9-57

25、6*x.7+432*x.5-120*x.3+9*x;elseifk=10w1=512*x.10-1280*x.8+1120*x.6-400*x.4+50*x.2-1;elseifk=11w1=1024*x.11-2816*x.9+2816*x.7-1232*x.5+220*x.3-11*x;elsew1=2048*x.12-6144*x.10+6912*x.8-3584*x.6+840*x.4-72*x.2+1;endfunctionw2=phi_yk(y,k)ifk=0w2=ones(size(y);elseifk=1w2=y;elseifk=2w2=2*y.2-1;elseifk=3w2=

26、4*y.3-3*y;elseifk=4w2=8*y.4-8*y.2+1;elseifk=5w2=16*y.5-20*y.3+5*y;elseifk=6w2=32*y.6-48*y.4+18*y.2-1;elseifk=7w2=64*y.7-112*y.5+56*y.3-7*y;elseifk=8w2=128*y.8-256*y.6+160*y.4-32*y.2+1;elseifk=9w2=256*y.9-576*y.7+432*y.5-120*y.3+9*y;elseifk=10w2=512*y.10-1280*y.8+1120*y.6-400*y.4+50*y.2-1;elseifk=11w

27、2=1024*y.11-2816*y.9+2816*y.7-1232*y.5+220*y.3-11*y;elsew2=2048*y.12-6144*y.10+6912*y.8-3584*y.6+840*y.4-72*y.2+1;endfunctiony=fun_phi1(x)%globali;globalj;y=phi_xk(x,i).*phi_xk(x,j);functionw=rho_phi(x,y)globali;globalj;w=obj(x,y).*phi_xk(x,i).*phi_yk(y,j);functiony=fun_phi(x)%globali;y=phi_xk(x,i).

28、*obj1(x);functionS=squar_approx(a,b,n)globali;globalj;ifnargin<3 n=1;endPhi2=zeros(n+1);fori=0:nfor j=0:n;Phi2(i+1,j+1)=quad(fun_phi1,a,b)-quad2d(rho_phi,a,b,a,(x)x); endendPhiF=zeros(n+1,1);fori=0:nPhiF(i+1)=quad(fun_phi,a,b);endS=Phi2PhiF;S=squar_approx(0,1,5);w=S'm,l=size(w);x = linspace(0

29、, 1, 5);U = 0*x;forj = 1:length(x)fork = 1:lU(j) = U(j) + w(k)*phi_xk(x(j),k-1);endendf1=exp(-x);aa=norm(U-f1)/norm(f1)fun= 'exp(-x)'fplot(fun,0,1)hold onplot(x,U,'o:')title( '真解与解析解的比较 ' );xlabel(' 变量 x' );ylabel(' 变量 x对应的值 ' );legend(' 真解','数值解 ');cc=U' f1' U'-f1'运行结果：aa =0.000725420168197738cc =1.0015012790446410.001501279044637370.9487402701690730.9487294800164371.07901526354981e-0050.8995614191959270.900087626252259-0.0005262070563323280.85342600033599