離散fredholm積分方程.doc

1 第一類Fredholm積分方程，具有形式如下：, (1)其中核函數(shù)和自由項(xiàng)為已知函數(shù)，是未知函數(shù)。此類積分方程雖然形式簡單，但其求解卻比較困難，所以這類方程在下文將做詳細(xì)介紹。2 第二類Fredholm積分方程，具有如下的形式：, (2)離散積分方程的數(shù)值方法有很多種，比如可以用復(fù)化梯形公式、復(fù)化辛普森公式等，這里我們利用復(fù)化梯形公式來進(jìn)行離散。一、復(fù)化梯形公式離散過程如下： 下面具體給出復(fù)化梯形公式對第二類積分方程的一般離散過程。 最后對變量進(jìn)行離散，將區(qū)間等分為份，步長為，同時(shí)忽略積分公式誤差項(xiàng)：其中 得到線性方程組其中，再對上述方程進(jìn)行數(shù)值求解，即可。例：求解積分方程，其解析解為代碼如下：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,各列分別表示真解、數(shù)值解、最小二乘解、正則解%a,b表示積分區(qū)間a,b%n表示將區(qū)間n等分%m表示正則參數(shù)的取值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維零向量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);%計(jì)算方程的真解endA(:,1)=A(:,1)/2;A(:,n+1)=A(:,n+1)/2;A=h*A;A=eye(n+1,n+1)-A;w2=Ag;%得到的數(shù)值解aa=norm(w1-w2)/norm(w1); %相對誤差bb=norm(w1-w2); %絕對誤差cc=w1 w2;plot(x,w1,b+)%真解hold onplot(x,w2,r*)%數(shù)值解%axis(0 1 -100 100);%設(shè)置坐標(biāo)軸title(數(shù)值解與真解的比較);%加圖形標(biāo)題xlabel(變量y);%加x軸說明ylabel(y對應(yīng)的解);%加y軸說明運(yùn)行結(jié)果： fww(0,1,50)aa = 0.000178436779942824 %相對誤差bb = 0.000693865370887685 %絕對誤差二、辛普森公式離散過程如下：下面給出復(fù)化梯形公式對第二類積分方程的一般離散過程。由于辛普森公式中取到中點(diǎn)的值，所以我們在區(qū)間上取個(gè)點(diǎn)。最后對變量進(jìn)行離散，將區(qū)間等分為份，步長為，同時(shí)忽略積分公式誤差項(xiàng)：其中 得到線性方程組其中，再對上述方程進(jìn)行數(shù)值求解，即可。例：求解積分方程，其解析解為代碼如下：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 quadrature 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(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,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 = 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*)%數(shù)值解%axis(0 1 -100 100);%設(shè)置坐標(biāo)軸title(數(shù)值解與真解的比較);%加圖形標(biāo)題xlabel(變量y);%加x軸說明ylabel(y對應(yīng)的解);%加y軸說明運(yùn)行結(jié)果：Enter number of pannels: 50aa = 6.9309212581323e-009aa = 2.69514294309828e-008三、高斯勒讓德離散過程如下：關(guān)于定積分，如果，則關(guān)于權(quán)函數(shù)正交多項(xiàng)式就是這時(shí)Gauss型積分公式的節(jié)點(diǎn)就取為上述多項(xiàng)式的零點(diǎn)，相應(yīng)的Gauss型積分公式為下面給出高斯公式對第二類積分方程的一般離散過程。在上Fredholm的積分公式為：第二類Fredholm積分方程可以化為：即高斯勒讓德型積分公式的積分區(qū)間為，而對于一般的區(qū)間上的積分需要作變量替換得到例：求積分方程其解析解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;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 %數(shù)值解plot(s,u);運(yùn)行結(jié)果：N=?5N = 5u = 0.999983504291541 1.00004348065192 0.999891437587449 1.00004348065192 0.999983504291541四、克倫肖柯蒂斯（Clenshaw-curtis）離散過程如下：五、高斯羅巴托（Gauss-Lobatto）離散： Gauss-Lobatto求積公式的表達(dá)式如下：Gauss-Lobatto求積公式的系數(shù)和余項(xiàng)分別為：其中，為的零點(diǎn)；為次Legendre(勒讓德)多項(xiàng)式六、伽遼金（Galerkin）法離散：設(shè)為空間內(nèi)的一個(gè)完備正交系，則當(dāng)充分大時(shí)，有其中為的逼近函數(shù)。將上式帶入得兩邊分別對求內(nèi)積，得即：得：其中例：求解第二類積分方程，其解析解為代碼如下：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-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.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; 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-2816*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).*obj1(x);function S=squar_approx(a,b,n) global i;global j; if nargin3 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 = 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對應(yīng)的值);legend(真解,數(shù)值解);cc=U f1 U-f1運(yùn)行結(jié)果：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.0004014826330763780.656714998388