fredholm,离散积分方程

1第一类Fredholm积分方程，具有形式如下：baxfdssysxk)()(),(,bxa(1)其中核函数),(sxK和自由项)(xf为已知函数，)(sy是未知函数。此类积分方程虽然形式简单，但其求解却比较困难，所以这类方程在下文将做详细介绍。2第二类Fredholm积分方程，具有如下的形式：baxfdssysxkxy)()(),()(,bxa(2)离散积分方程的数值方法有很多种，比如可以用复化梯形公式、复化辛普森公式等，这里我们利用复化梯形公式来进行离散。一、复化梯形公式离散过程如下：)]()(2)([2)(1bfxfafhdxxfnkkba下面具体给出复化梯形公式对第二类积分方程的一般离散过程。],[)()(),(12)()](),(21)(),()(),()(),(21[)(),(12)()](),()(),([2)](),()(),([2)](),()(),([2)(),()(),()(),()(),()(),(21100211111100112110baxgfxkabhsysxksysxksysxksysxkhfxkabhsysxksysxkhsysxksysxkhsysxksysxkhdssysxkdssysxkdssysxkdssysxkdssysxknniinnnniiiissssssssbannii－最后对变量x进行离散，将区间],[ba等分为n份，步长为nabh，同时忽略积分公式误差项：)](),(21)(),()(),()(),(21[)()(1100nniiiiiiiisysxksysxksysxksysxkhxyxg其中ni,2,1,0得到线性方程组nngAf其中)](),(),(),([210nnsysysysyf，)](,),(),(),([210nnxgxgxgxgg),(211),(),(21),(21),(1),(21),(21),(),(211101110101000nnnnnnyxhkyxhkyxkhyxhkyxhkyxhkyxhkyxhkyxhkA再对上述方程进行数值求解，即可。例：求解积分方程10223)1(85254)()11()(xxxxdssyxsxy，其解析解为2)1()(xxy代码如下：functionK=K(x,y)K=1/(1+y)-x;functionw1=fun1(x)w1=1./((1+x).*(1+x));functionf=f(x)f=(4*x.*x.*x+5*x.*x-2*x+5)./(8*(x+1).*(x+1));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=1:n+1forj=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=A\g;%得到的数值解aa=norm(w1-w2)/norm(w1);%相对误差bb=norm(w1-w2);%绝对误差cc=[w1w2];plot(x,w1,'b+')%真解holdonplot(x,w2,'r*')%数值解%axis([01-100100]);%设置坐标轴title('数值解与真解的比较');%加图形标题xlabel('变量y');%加x轴说明ylabel('y对应的解');%加y轴说明运行结果：fww(0,1,50)aa=0.000178436779942824%相对误差bb=0.000693865370887685%绝对误差00.10.20.30.40.50.60.70.80.910.20.30.40.50.60.70.80.911.11.2数值解与真解的比较变量yy对应的解二、辛普森公式离散过程如下：)]()2(4)([6)(bfbafafhdxxfba下面给出复化梯形公式对第二类积分方程的一般离散过程。由于辛普森公式中取到中点的值，所以我们在区间],[ba上取12n个点。],[)()(2880)()](),()(),(4)(),(2)(),(4)(),([6)(2880)()](),()2()2,(4)(),([6)](),()2()2,(4)(),([6)](),()2()2,(4)(),([6)(),()(),()(),()(),()4(4221212221100)4(41111222121211110100012110baxgfhabsysxksysxksysxksysxksysxkhfhabsysxkssyssxksysxkhsysxkssyssxksysxkhsysxkssyssxksysxkhdssysxkdssysxkdssysxkdssysxknnnnkkkkkkkkssssssbakk最后对变量x进行离散，将区间],[ba等分为n2份，步长为nabh2，同时忽略积分公式误差项：)](),()(),(4)(),(2)(),(4)(),([6)()(221212221100nninniiiiiisysxksysxksysxksysxksysxkhxyxg其中ni2,,2,1,0得到线性方程组nngAf22其中)](,),(),(),([22102nnsysysysyf，)](,),(),(),([22102nnxgxgxgxgg),(611),(31),(32),(61),(61),(31),(321),(61),(61),(31),(32),(611),(),(2),(4),(),(),(2),(4),(),(),(2),(4),(6222212022121110120201000222212022121110120201000nnnnnnnnnnnnnnyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxkyxkyxkyxkyxkyxkyxkyxkyxkyxkyxkyxkhIA再对上述方程进行数值求解，即可。例：求解积分方程10223)1(85254)()11()(xxxxdssyxsxy，其解析解为2)1()(xxy代码如下：functionv=knl(x,t)v=1/(1+t)-x;functionf=fnc(x)f=(4*x.*x.*x+5*x.*x-2*x+5)./(8*(x+1).*(x+1));functiony=inteqn(t,kernel,fun,coef)%%Inputs%tevaluationpointsofthequadraturerule%kernelkernelfunctionK%funfunctionf%coefquadraturerulecoefficients%Output%ydiscretesolutionvaluesattn=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);forj=1:nA(:,j)=-coef(j)*K(:,j);A(j,j)=1.0+A(j,j);endy=A\f;k=input('Enternumberofpannels:');x=linspace(0,1,k+1);x=x';%evenlyspacedknots%Note:thisxcanbereplacedbyanypartitionof[0,1]y=zeros(length(x),1);%discreteapproximationatx%useSimpson’sn=2*k+1;%numberofpointsinSimpson’srulecoef=zeros(n,1);%coeficientsinSimpson’trulet=zeros(n,1);%pointsinSimpson’srule%generateSimpson’srulecoefficientsandevaluationpointsfori=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);%solvetheintegralequationyt=inteqn(t,@knl,@fnc,coef);%discreteapproximationtoy(x)atthepartitionxy=yt(1:2:n);%checkresultsyexact=1./((1+x).*(1+x));aa=norm(y-yexact)/norm(y)%相对误差aa=norm(y-yexact)%绝对误差cc=[yyexact]plot(x,y,'b+')%真解holdonplot(x,yexact,'r*')%数值解%axis([01-100100]);%设置坐标轴title('数值解与真解的比较');%加图形标题xlabel('变量y');%加x轴说明ylabel('y对应的解');%加y轴说明运行结果：Enternumberofpannels:50aa=6.9309212581323e-009aa=2.69514294309828e-00800.10.20.30.40.50.60.70.80.910.20.30.40.50.60.70.80.911.11.2数值解与真解的比较变量yy对应的解三、高斯－勒让德离散过程如下：关于定积分11-)()(dxxfx，如果1)(],1,1[],[xba，则关于权函数1)(x正交多项式就是nnmnndxxdnxL)1(!21)(2这时Gauss型积分公式的节点就取为上述多项式的零点，相应的Gauss型积分公式为nkkkxfAdxxf111-)()(下面给出高斯公式对第二类积分方程的一般离散过程。在]1,1-[上Fredholm的积分公式为：)()(),()(),(111-xgsysxkAdssysxknjjjj第二类Fredholm积分方程可以化为：njjjijiisysxkAxyxg1)(),()()(即)()()()()()(),(1),(),(),(),(1),(),(),(),(12121221122221211212111nnnnnnnnnnnxgxgxgsysysysxkAsxkAsxkAsxkAsxkAsxkAsxkAsxkAsxkA高斯－勒让德型积分公式的积分区间为]1,1[，而对于一般的区间],[ba上的积分dxxfba)(需要作变量替换t