fredholm,离散积分方程

1、1 第一类Fredholm积分方程，具有形式如下：, (1)其中核函数和自由项为已知函数，是未知函数。此类积分方程虽然形式简单，但其求解却比较困难，所以这类方程在下文将做详细介绍。2 第二类Fredholm积分方程，具有如下的形式：, (2)离散积分方程的数值方法有很多种，比如可以用复化梯形公式、复化辛普森公式等，这里我们利用复化梯形公式来进行离散。一、复化梯形公式离散过程如下： 下面具体给出复化梯形公式对第二类积分方程的一般离散过程。 最后对变量进行离散，将区间等分为份，步长为，同时忽略积分公式误差项：其中 得到线性方程组其中，再对上述方程进行数值求解，即可。例：求解积分方程，其解析解为代码

2、如下：function K=K(x,y)K = 1/(1+y) - x;function w1=fun1(x)w1=1./(1+x).*(1+x);function f=f(x)f = (4*x.*x.*x + 5*x.*x - 2*x + 5)./(8*(x+1).*(x+1);function w5=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);%初始化

3、矩阵A为n+1阶零矩阵g=zeros(n+1,1);%初始化列向量g为n+1维零向量w1=zeros(n+1,1);%初始化列向量w1为n+1维零向量for i=1:n+1 for j=1:n+1 A(i,j)=K(x(i),y(j); end g(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

4、(x,w1,'b+')%真解hold onplot(x,w2,'r*')%数值解%axis(0 1 -100 100);%设置坐标轴title('数值解与真解的比较');%加图形标题xlabel('变量y');%加x轴说明ylabel('y对应的解');%加y轴说明运行结果：>> fww(0,1,50)aa = 0.000178436779942824 %相对误差bb = 0.000693865370887685 %绝对误差二、辛普森公式离散过程如下：下面给出复化梯形公式对第二类积分方程的一般离散过程。

5、由于辛普森公式中取到中点的值，所以我们在区间上取个点。最后对变量进行离散，将区间等分为份，步长为，同时忽略积分公式误差项：其中 得到线性方程组其中，再对上述方程进行数值求解，即可。例：求解积分方程，其解析解为代码如下：function v = knl(x,t)v = 1/(1+t) - x;function f = fnc(x)f = (4*x.*x.*x + 5*x.*x - 2*x + 5)./(8*(x+1).*(x+1);function y = inteqn(t, kernel, fun, coef)% Inputs% t evaluation points of the quadr

6、ature rule% kernel kernel function K% fun function f% coef quadrature rule coefficients% Output% y discrete solution values at tn = length(t);f = feval(fun, t);%for j=1:nfor i = 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(

7、j,j);endy = Af;k = input('Enter number of pannels: ');x = linspace(0,1,k+1); x = x' % evenly spaced knots% Note: this x can be replaced by any partition of 0,1y = zeros(length(x),1); % discrete approximation at x% use Simpsonsn = 2*k + 1; % number of points in Simpsons rulecoef = zeros(n

8、,1); % coeficients in Simpsont rulet = zeros(n,1); % points in Simpsons rule% generate Simpsons rule coefficients and evaluation pointsfor i=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

9、 = inteqn(t, knl, fnc, coef);% discrete approximation to y(x) at the partition xy = yt(1:2:n);% check resultsyexact = 1./(1+x).*(1+x); aa=norm(y-yexact)/norm(y) %相对误差aa=norm(y - yexact) %绝对误差cc=y yexactplot(x,y,'b+')%真解hold onplot(x,yexact,'r*')%数值解%axis(0 1 -100 100);%设置坐标轴title(

10、9;数值解与真解的比较');%加图形标题xlabel('变量y');%加x轴说明ylabel('y对应的解');%加y轴说明运行结果：Enter number of pannels: 50aa = 6.9309212581323e-009aa = 2.69514294309828e-008三、高斯勒让德离散过程如下：关于定积分，如果，则关于权函数正交多项式就是这时Gauss型积分公式的节点就取为上述多项式的零点，相应的Gauss型积分公式为下面给出高斯公式对第二类积分方程的一般离散过程。在上Fredholm的积分公式为：第二类Fredholm积分方程可以

11、化为：即高斯勒让德型积分公式的积分区间为，而对于一般的区间上的积分需要作变量替换得到例：求积分方程其解析解function x,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) ;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

12、'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.999983504291541 1.00004348065192 0.999891437587449 1.00004348065192 0.999983504291541四、克伦肖柯蒂斯（Clenshaw-curtis）离散过程如下：五、高斯罗巴托（Gauss-Lobatto）离散： Gauss-Lobatto求积公式的表达式如下：G

13、auss-Lobatto求积公式的系数和余项分别为：其中，为的零点；为次Legendre(勒让德)多项式六、伽辽金（Galerkin）法离散：设为空间内的一个完备正交系，则当充分大时，有其中为的逼近函数。将上式带入得两边分别对求内积，得即：得：其中例：求解第二类积分方程，其解析解为代码如下：function w=obj(x,y)w=exp(-x-y);function w=obj1(x)w=(exp(-x)+exp(-3*x)/2;function w1=phi_xk(x,k)if k=0 w1=ones(size(x);elseif k=1 w1=x; elseif k=2 w1=2*x.2

14、-1; elseif k=3 w1=4*x.3-3*x; elseif k=4 w1=8*x.4-8*x.2+1; elseif k=5 w1=16*x.5-20*x.3+5*x; elseif k=6 w1=32*x.6-48*x.4+18*x.2-1; elseif k=7 w1=64*x.7-112*x.5+56*x.3-7*x; elseif k=8 w1=128*x.8-256*x.6+160*x.4-32*x.2+1; elseif k=9 w1=256*x.9-576*x.7+432*x.5-120*x.3+9*x; elseif k=10 w1=512*x.10-1280*x.

15、8+1120*x.6-400*x.4+50*x.2-1; elseif k=11 w1=1024*x.11-2816*x.9+2816*x.7-1232*x.5+220*x.3-11*x; else w1=2048*x.12-6144*x.10+6912*x.8-3584*x.6+840*x.4-72*x.2+1;endfunction w2=phi_yk(y,k)if k=0 w2=ones(size(y);elseif k=1 w2=y; elseif k=2 w2=2*y.2-1; elseif k=3 w2=4*y.3-3*y; elseif k=4 w2=8*y.4-8*y.2+1;

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

17、6*y.9+2816*y.7-1232*y.5+220*y.3-11*y; else w2=2048*y.12-6144*y.10+6912*y.8-3584*y.6+840*y.4-72*y.2+1;endfunction y=fun_phi1(x) %global i;global j;y=phi_xk(x,i).*phi_xk(x,j);function w=rho_phi(x,y)global i;global j;w=obj(x,y).*phi_xk(x,i).*phi_yk(y,j);function y=fun_phi(x) %global i;y=phi_xk(x,i).*ob

18、j1(x);function S=squar_approx(a,b,n) global i;global j; if nargin<3 n=1;endPhi2=zeros(n+1); for i=0:n for j=0:n; Phi2(i+1,j+1)=quad(fun_phi1,a,b)-quad2d(rho_phi,a,b,a,(x)x); end endPhiF=zeros(n+1,1); for i=0:n PhiF(i+1)=quad(fun_phi,a,b);end S=Phi2PhiF;S=squar_approx(0,1,5);w=S'm,l=size(w);x

19、= linspace(0, 1, 5);U = 0*x;for j = 1:length(x)for k = 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('真解','数值解

20、');cc=U' f1' U'-f1'运行结果：aa = 0.000725420168197738cc =1.00150127904464 1 0.001501279044637370.948740270169073 0.948729480016437 1.07901526354981e-0050.899561419195927 0.900087626252259 -0.0005262070563323280.853426000335993 0.853939665623535 -0.0005136652875424860.809908936106422 0.810157734932427 -0.000248798826004260.768684772899146 0.768620526593736 6.42463054101317e-0050.729513656549289 0.729212952525235 0.0003007040240542440.692227307903594 0.691825825270517 0.0004014