2022年fredholm,離散積分方程

1、1 第一類 fredholm 積分方程，具有形式如下：baxfdssysxk)()(),(,bxa (1) 其中核函數(shù)),(sxk和自由項)(xf為已知函數(shù)，)(sy是未知函數(shù)。此類積分方程雖然形式簡單， 但其求解卻比較困難， 所以這類方程在下文將做詳細(xì)介紹。2 第二類 fredholm 積分方程，具有如下的形式：baxfdssysxkxy)()(),()(,bxa (2) 離散積分方程的數(shù)值方法有很多種，比如可以用復(fù)化梯形公式、復(fù)化辛普森公式等，這里我們利用復(fù)化梯形公式來進(jìn)行離散。一、復(fù)化梯形公式離散過程如下：)()(2)(2)(1bfxfafhdxxfnkkba下面具體給出復(fù)化梯形公式對第

2、二類積分方程的一般離散過程。,)()(),(12)()(),(21)(),()(),()(),(21)(),(12)()(),()(),(2)(),()(),(2)(),()(),(2)(),()(),()(),()(),()(),(21100211111100112110baxgfxkabhsysxksysxksysxksysxkhfxkabhsysxksysxkhsysxksysxkhsysxksysxkhdssysxkdssysxkdssysxkdssysxkdssysxknniinnnniiiissssssssbannii最后對變量x進(jìn)行離散，將區(qū)間,ba等分為n份， 步長為nabh，

3、同時忽略積分公式誤差項：精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 1 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 1 頁，共 14 頁 - - - - - - - - -)(),(21)(),()(),()(),(21)()(1100nniiiiiiiisysxksysxksysxksysxkhxyxg其中ni,2,1 ,0得到線性方程組nngaf其中)(),(),(),(210nnsysysysyf，)(,),(),(),(210nnxgxgx

4、gxgg),(211),(),(21),(21),(1),(21),(21),(),(211101110101000nnnnnnyxhkyxhkyxkhyxhkyxhkyxhkyxhkyxhkyxhka再對上述方程進(jìn)行數(shù)值求解，即可。例：求解積分方程10223) 1(85254)()11()(xxxxdssyxsxy， 其解析解為2)1()(xxy代碼如下：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)

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+1for j=1:n+1 a(i,j)=k(x(i),

6、y(j);end精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 2 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 2 頁，共 14 頁 - - - - - - - - - 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; % 得到的數(shù)值解aa=norm(w1-w2)/norm(w1); %

7、相對誤差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軸說明運行結(jié)果： fww(0,1,50) aa = 0.000178436779942824% 相對誤差bb = 0.000693865370887685 % 絕對誤差00.10.20.30.40.50.60.70.80.910.20.30.40.5

8、0.60.70.80.911.11.2數(shù) 值 解 與 真 解 的 比 較變 量 yy對應(yīng)的解精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 3 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 3 頁，共 14 頁 - - - - - - - - -二、辛普森公式離散過程如下：)()2(4)(6)(bfbafafhdxxfba下面給出復(fù)化梯形公式對第二類積分方程的一般離散過程。由于辛普森公式中取到中點的值，所以我們在區(qū)間,ba上取12n個點。,)()(288

9、0)()(),()(),(4)(),(2)(),(4)(),(6)(2880)()(),()2()2,(4)(),(6)(),()2()2,(4)(),(6)(),()2()2,(4)(),(6)(),()(),()(),()(),()4(4221212221100)4(41111222121211110100012110baxgfhabsysxksysxksysxksysxksysxkhfhabsysxkssyssxksysxkhsysxkssyssxksysxkhsysxkssyssxksysxkhdssysxkdssysxkdssysxkdssysxknnnnkkkkkkkkssssss

10、bakk最 后對變量x進(jìn) 行離 散，將區(qū)間,ba等分 為n2份，步 長為nabh2，同時忽略積分公式誤差項：)(),()(),(4)(),(2)(),(4)(),(6)()(221212221100nninniiiiiisysxksysxksysxksysxksysxkhxyxg其中ni2,2 ,1 ,0得到線性方程組nngaf22其中)(,),(),(),(22102nnsysysysyf，)(,),(),(),(22102nnxgxgxgxgg精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 4 頁，共 14 頁 - - - - - - - - -

11、精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 4 頁，共 14 頁 - - - - - - - - -),(611),(31),(32),(61),(61),(31),(321),(61),(61),(31),(32),(611),(),(2),(4),(),(),(2),(4),(),(),(2),(4),(6222212022121110120201000222212022121110120201000nnnnnnnnnnnnnnyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxhkyxkyxkyxky

12、xkyxkyxkyxkyxkyxkyxkyxkyxkhia再對上述方程進(jìn)行數(shù)值求解，即可。例：求解積分方程10223) 1(85254)()11()(xxxxdssyxsxy，其解析解為2)1 ()(xxy代碼如下：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 quadra

13、ture 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:n精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - -

14、 第 5 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 5 頁，共 14 頁 - - - - - - - - -a(:,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

15、 0,1y = zeros(length(x),1); % discrete approximation at x% use simpson sn = 2*k + 1; % number of points in simpson s rulecoef = zeros(n,1); % coeficients in simpson t rulet = zeros(n,1); % points in simpson s rule% generate simpson s rule coefficients and evaluation pointsfor i=1:kcoef(2*i-1:2*i+1)

16、= 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) %

17、相對誤差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軸說明運行結(jié)果：enter number of pannels: 50 aa = 精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 6 頁，共 14 頁 - - - - - - - -

18、-精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 6 頁，共 14 頁 - - - - - - - - -6.9309212581323e-009 aa = 2.69514294309828e-008 00.10.20.30.40.50.60.70.80.910.20.30.40.50.60.70.80.911.11.2數(shù) 值 解 與 真 解 的 比 較變 量 yy對應(yīng)的解三、高斯勒讓德離散過程如下：關(guān)于定積分11 -)()(dxxfx，如果1)(,1 , 1,xba，則關(guān)于權(quán)函數(shù)1)(x正交多項式就是nnmnndxxdnxl)1(!21)(2這時

19、gauss型積分公式的節(jié)點就取為上述多項式的零點，相應(yīng)的gauss型積分公式為nkkkxfadxxf111-)()(下面給出高斯公式對第二類積分方程的一般離散過程。精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 7 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 7 頁，共 14 頁 - - - - - - - - -在 1 , 1- 上 fredholm 的積分公式為：)()(),()(),(111-xgsysxkadssysxknjjjj第二類 fre

20、dholm 積分方程可以化為：njjjijiisysxkaxyxg1)(),()()(即)()()()()()(),(1),(),(),(),(1),(),(),(),(12121221122221211212111nnnnnnnnnnnxgxgxgsysysysxkasxkasxkasxkasxkasxkasxkasxkasxka高斯勒讓德型積分公式的積分區(qū)間為1 , 1， 而對于一般的區(qū)間,ba上的積分dxxfba)(需要作變量替換tababx22得到11)22(2)(dttababfabdxxfba例：求積分方程)1arctan()1arctan()()(111)(112sxdssfsx

21、xy其解析解1)(xyfunction 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

22、-s)+atan(1+s);ww=cond(k*w)精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 8 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 8 頁，共 14 頁 - - - - - - - - -u=ag % 數(shù)值解plot(s,u);運行結(jié)果：n=?5 n = 5 u = 0.999983504291541 1.00004348065192 0.999891437587449 1.00004348065192 0.99998350429154

23、1 四、克倫肖柯蒂斯（ clenshaw-curtis）離散過程如下：五、高斯羅巴托（ gauss-lobatto ）離散：gauss-lobatto 求積公式的表達(dá)式如下：)()1()1() 1(2)(1211-frxfaffnndxxfnnjjjgauss-lobatto 求積公式的系數(shù)和余項分別為：2)()1(2jnjxpnna) 1 , 1(, )()!2)(12()!1(2) 1()2(34123nnnfnnnnnfr其中，jx為)(xpn的零點；)(xpn為n次legendre(勒讓德)多項式六、伽遼金（ galerkin ）法離散：設(shè))(xi為,2bal空間內(nèi)的一個完備正交系，,

24、)(2balxk， 則當(dāng)n精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 9 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 9 頁，共 14 頁 - - - - - - - - -充分大時，有)()(0 xkxiniin其中)(xn為)(x的逼近函數(shù)。將上式帶入baxgdssysxkxy)()(),()(得bainiiiniixgdsssxkkxk)()(),()(00兩邊分別對)(xj求內(nèi)積，得)(),()(, )(),()(0 xxgxssxkxkj

25、jbaibainii即：dxxxgkdsdxxssxkdxxxbajniibababajiji)()()()(),()()(0得：nnijbbbbkkkka210210其中ijabababajijidsdxxssxkdxxx)()(),()()(dxxxgbbajj)()(例：求解第二類積分方程10, )(21)()(03xxxyxxeedssyexy，其解析解為xexy)(代碼如下：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精品學(xué)習(xí)資

26、料 可選擇p d f - - - - - - - - - - - - - - 第 10 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 10 頁，共 14 頁 - - - - - - - - - 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

27、*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-

28、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

29、k=8精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 11 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 11 頁，共 14 頁 - - - - - - - - - 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;e

30、lseif 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) %

31、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: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); for i=0:n phif(i+1)=quad(fun_phi,a,b);ends=phi2phif;s=squar_approx(0,1,5);w=s;

32、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);endend精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 12 頁，共 14 頁 - - - - - - - - -精品學(xué)習(xí)資料 可選擇p d f - - - - - - - - - - - - - - 第 12 頁，共 14 頁 - - - - - - - - -f1=exp(-x);aa=norm(u-f1)/norm(f1)fun

33、= 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 運行結(jié)果：aa = 0.000725420168197738 cc = 1.00150127904464 1 0.00150127904463737 0.948740270169073 0.948729480016437 1.07901526354981e-005 0.899561419195927 0.900087626252259 -0.000526207056332328 0.853426000335993 0.853939665623535 -0.000513665287542486 0.809908936106422 0.810157734932427 -0.00024879882600426 0.768684772899146 0.768620526593736 6.42463054101317e-005 0.729513656549289 0.729212952525235 0.000300704024054244 0.6922