matlab代码--数值积分.doc

数值积分1. CombineTraprl复合梯形公式求积分function I,step = CombineTraprl(f,a,b,eps)%f 被积函数%a,b 积分上下限%eps 精度%I 积分结果%step 积分的子区间数if(nargin =3) eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h;while abs(I2-I1)eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f),x)+subs(sym(f),findsym(sym(f),x1); endendI=I2;step=n; 2. IntSimpson用辛普森系列公式求积分function I,step = IntSimpson(f,a,b,type,eps)%type = 1 辛普森公式%type = 2 辛普森3/8公式%type = 3 复合辛普森公式if(type=3 & nargin=4) eps=1.0e-4; %缺省精度为0.0001endI=0;switch type case 1, I=(b-a)/6)*(subs(sym(f),findsym(sym(f),a)+. 4*subs(sym(f),findsym(sym(f),(a+b)/2)+. subs(sym(f),findsym(sym(f),b); step=1; case 2, I=(b-a)/8)*(subs(sym(f),findsym(sym(f),a)+. 3*subs(sym(f),findsym(sym(f),(2*a+b)/3)+ . 3*subs(sym(f),findsym(sym(f),(a+2*b)/3)+subs(sym(f),findsym(sym(f),b); step=1; case 3, n=2; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h; while abs(I2-I1)eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/6)*(subs(sym(f),findsym(sym(f),x)+. 4*subs(sym(f),findsym(sym(f),(x+x1)/2)+. subs(sym(f),findsym(sym(f),x1); end end I=I2; step=n;end 3. NewtonCotes用牛顿科茨系列公式求积分function I = NewtonCotes(f,a,b,type)%type = 1 科茨公式%type = 2 牛顿科茨六点公式%type = 3 牛顿科茨七点公式I=0;switch type case 1, I=(b-a)/90)*(7*subs(sym(f),findsym(sym(f),a)+. 32*subs(sym(f),findsym(sym(f),(3*a+b)/4)+. 12*subs(sym(f),findsym(sym(f),(a+b)/2)+. 32*subs(sym(f),findsym(sym(f),(a+3*b)/4)+7*subs(sym(f),findsym(sym(f),b); case 2, I=(b-a)/288)*(19*subs(sym(f),findsym(sym(f),a)+. 75*subs(sym(f),findsym(sym(f),(4*a+b)/5)+. 50*subs(sym(f),findsym(sym(f),(3*a+2*b)/5)+. 50*subs(sym(f),findsym(sym(f),(2*a+3*b)/5)+. 75*subs(sym(f),findsym(sym(f),(a+4*b)/5)+19*subs(sym(f),findsym(sym(f),b); case 3, I=(b-a)/840)*(41*subs(sym(f),findsym(sym(f),a)+. 216*subs(sym(f),findsym(sym(f),(5*a+b)/6)+. 27*subs(sym(f),findsym(sym(f),(2*a+b)/3)+. 272*subs(sym(f),findsym(sym(f),(a+b)/2)+. 27*subs(sym(f),findsym(sym(f),(a+2*b)/3)+. 216*subs(sym(f),findsym(sym(f),(a+5*b)/6)+41*subs(sym(f),findsym(sym(f),b);end4. IntGauss用高斯公式求积分function I = IntGauss(f,a,b,n,AK,XK)if(n5 & nargin = 4) AK = 0; XK = 0;else XK1=(b-a)/2)*XK+(a+b)/2); I=(b-a)/2)*sum(AK.*subs(sym(f),findsym(f),XK1);endta = (b-a)/2;tb = (a+b)/2;switch n case 0, I=2*ta*subs(sym(f),findsym(sym(f),tb); case 1, I=ta*(subs(sym(f),findsym(sym(f),ta*0.5773503+tb)+. subs(sym(f),findsym(sym(f),-ta*0.5773503+tb); case 2, I=ta*(0.55555556*subs(sym(f),findsym(sym(f),ta*0.7745967+tb)+. 0.55555556*subs(sym(f),findsym(sym(f),-ta*0.7745967+tb)+. 0.88888889*subs(sym(f),findsym(sym(f),tb); case 3, I=ta*(0.3478548*subs(sym(f),findsym(sym(f),ta*0.8611363+tb)+. 0.3478548*subs(sym(f),findsym(sym(f),-ta*0.8611363+tb)+. 0.6521452*subs(sym(f),findsym(sym(f),ta*0.3398810+tb). +0.6521452*subs(sym(f),findsym(sym(f),-ta*0.3398810+tb); case 4, I=ta*(0.2369269*subs(sym(f),findsym(sym(f),ta*0.9061793+tb)+. 0.2369269*subs(sym(f),findsym(sym(f),-ta*0.9061793+tb)+. 0.4786287*subs(sym(f),findsym(sym(f),ta*0.5384693+tb). +0.4786287*subs(sym(f),findsym(sym(f),-ta*0.5384693+tb)+. 0.5688889*subs(sym(f),findsym(sym(f),tb);end 5. IntGaussLada用高斯拉道公式求积分function I = IntGaussLada(f,a,b,n,AK,XK)if(n6 & nargin = 4) AK = 0; XK = 0;else XK1=(b-a)/2)*XK+(a+b)/2); I=(b-a)/2)*(2/n/n)*subs(sym(f),findsym(sym(f),a)+. sum(AK.*subs(sym(f),findsym(sym(f),XK1);endta = (b-a)/2;tb = (a+b)/2;switch n case 2, I=ta*0.5*subs(sym(f),findsym(sym(f),ta*(-1)+tb)+. 1.5*ta*subs(sym(f),findsym(sym(f),ta*(1/3)+tb); case 3, I=ta*(2/9)*subs(sym(f),findsym(sym(f),ta*(-1)+tb)+. 0.752806*subs(sym(f),findsym(sym(f),ta*(-0.289898)+tb)+. 1.024972*subs(sym(f),findsym(sym(f),ta*0.689898+tb); case 4, I=ta*(0.125*subs(sym(f),findsym(sym(f),ta*(-1)+tb)+. 0.657689*subs(sym(f),findsym(sym(f),-ta*0.575319+tb)+. 0.776387*subs(sym(f),findsym(sym(f),ta*0.181066+tb)+. 0.440925*subs(sym(f),findsym(sym(f),ta*0.822824+tb); case 5, I=ta*(0.08*subs(sym(f),findsym(sym(f),ta*(-1)+tb)+. 0.446207*subs(sym(f),findsym(sym(f),-ta*0.720480+tb)+. 0.623653*subs(sym(f),findsym(sym(f),-ta*0.167181+tb)+. 0.562712*subs(sym(f),findsym(sym(f),ta*0.446314+tb)+. 0.287427*subs(sym(f),findsym(sym(f),ta*0.885792+tb); end6. IntGaussLobato用高斯洛巴托公式求积分function I = IntGaussLobato(f,a,b,n,AK,XK)if(n7 & nargin = 4) AK = 0; XK = 0;else XK1=(b-a)/2)*XK+(a+b)/2); I=(b-a)/2)*(2/n/(n-1)*(subs(sym(f),findsym(sym(f),a)+. subs(sym(f),findsym(sym(f),b)+sum(AK.*subs(sym(f),findsym(sym(f),XK1);endta = (b-a)/2;tb = (a+b)/2;switch n case 3, I=ta*(1/3)*(subs(sym(f),findsym(sym(f),a)+. subs(sym(f),findsym(sym(f),b)+1.333333*subs(sym(f),findsym(sym(f),tb); case 4, I=ta*(1/6)*(subs(sym(f),findsym(sym(f),a)+. subs(sym(f),findsym(sym(f),b)+0.833333*. (subs(sym(f),findsym(sym(f),ta*0.447214+tb)+. subs(sym(f),findsym(sym(f),-ta*0.447214+tb); case 5, I=ta*(1/10)*(subs(sym(f),findsym(sym(f),a)+. subs(sym(f),findsym(sym(f),b)+0.544444*. (subs(sym(f),findsym(sym(f),ta*0.654654+tb)+. subs(sym(f),findsym(sym(f),-ta*0.654654+tb)+. 0.711111*subs(sym(f),findsym(sym(f),tb); case 6, I=ta*(1/15)*(subs(sym(f),findsym(sym(f),a)+. subs(sym(f),findsym(sym(f),b)+0.554858*. (subs(sym(f),findsym(sym(f),ta*0.285232+tb)+. subs(sym(f),findsym(sym(f),-ta*0.285232+tb)+. 0.378475*(subs(sym(f),findsym(sym(f),ta*0.765055+tb)+. subs(sym(f),findsym(sym(f),-ta*0.765055+tb); end7. IntSample用三次样条插值求积分function q = IntSample(func,a,b,n)format long;node_num = n + 3; %加上延拓的2个节点共n+3个节点h = (b-a)/n;for i=1:node_num y(i) = subs(sym(func), findsym(sym(func), a+i*h-2*h);endy_1 = (-3*y(1)+4*y(2)-y(3)/(2*h);y_N = (3*y(node_num)-4*y(node_num-1)+3*y(node_num-2)/(2*h);c = SubBSample(h,node_num-1,y,y_1,y_N); %用样条逼近法求出样条系数q = h*(c(1)+c(node_num+2)/24+h*(c(2)+c(node_num+1)/2+23*h*(c(3)+c(node_num)/24;for j=4:node_num-1 q = q+h*c(j);endformat short;function c = SubBSample(h,n,y,y_1,y_N)f0 = 0.0;c = zeros(n+3,1);b = zeros(n+1,1);A = diag(4*ones(n+1,1);I = eye(n+1,n+1);AL = I(2:n+1,:);zeros(1,n+1);AU = zeros(1,n+1);I(1:n,:);A = A+AL+AU; %形成系数矩阵for i=2:n b(i,1) = 6*y(i);endb(1) = 6*y(1)+2*h*y_1;b(n+1) = 6*y(n+1)-2*h*y_N;d = followup(A,b); %用追赶法求出系数c(2:n+2) = d;c(1) = c(2) - 2*h*y_1; %c(-1)c(n+3) = c(3)+2*h*y_N; %c(n+1)8. IntPWC用抛物插值求积分function q = IntPWC(X,Y,n)format long;h = zeros(n-1,1);lamda = zeros(n,1);L = zeros(n-2,1);R = zeros(n-2,1);Dlta = zeros(n-2,1);h(1) = X(2)-X(1);for j=1:n-2 %计算每个子区间的积分参数 h(j+1) = X(j+2)-X(j+1); lamda(j) = h(j)/h(j+1); Dlta(j) = h(j+1)/h(j); L(j) = (Dlta(j)*Dlta(j)*Y(j)/(1+Dlta(j)-Dlta(j)*Y(j+1)+. +Dlta(j)*Y(j+2)/(1+Dlta(j); R(j) = (lamda(j)*lamda(j)*Y(j+2)/(1+lamda(j)-lamda(j)*Y(j+1)+. +lamda(j)*Y(j)/(1+lamda(j);endq = h(1)*(3*Y(1)+3*Y(2)-R(1)/6+h(n-1)*(3*Y(n-1)+3*Y(n)-L(n-2)/6;for k=2:n-2 q = q + h(k)*(3*Y(k)+3*Y(k+1)-0.5*R(k)-0.5*L(k-1)/6;endformat short;9. IntGaussLager用高斯-拉盖尔公式求积分function I = IntGaussLager(f,n,AK,XK)if(n6 & nargin = 2) AK = 0; XK = 0;else I=sum(AK.*subs(sym(f),findsym(sym(f),XK);endswitch n case 2, I=0.853553*subs(sym(f),findsym(sym(f),-0.585786)+. 0.146447*subs(sym(f),findsym(sym(f),3.414214); case 3, I=0.711093*subs(sym(f),findsym(sym(f),0.415575)+. 0.278518*subs(sym(f),findsym(sym(f),2.294280)+. 0.0103893*subs(sym(f),findsym(sym(f),6.289945); case 4, I=0.603154*subs(sym(f),findsym(sym(f),0.322548)+. 0.357419*subs(sym(f),findsym(sym(f),1.745761)+. 0.0388879*subs(sym(f),findsym(sym(f),4.536620)+. 0.000539295*subs(sym(f),findsym(sym(f),9.395071); case 5, I=0.521756*subs(sym(f),findsym(sym(f),0.263560)+. 0.398667*subs(sym(f),findsym(sym(f),1.413403)+. 0.0759424*subs(sym(f),findsym(sym(f),3.596426)+. 0.00361176*subs(sym(f),findsym(sym(f),7.085810)+. 0.0000233700*subs(sym(f),findsym(sym(f),12.640801); end10. IntGaussHermite用高斯-埃尔米特公式求积分function I = IntGaussHermite(f,n,AK,XK)if(n6 & nargin = 2) AK = 0; XK = 0;else I=sum(AK.*subs(sym(f),findsym(sym(f),XK);endswitch n case 2, I=0.886227*(subs(sym(f),findsym(sym(f),-0.707107)+. subs(sym(f),findsym(sym(f),0.707107); case 3, I=1.181636*subs(sym(f),findsym(sym(f),0)+. 0.295409*(subs(sym(f),findsym(sym(f),1.224745)+. subs(sym(f),findsym(sym(f),-1.224745); case 4, I=0.544444*(subs(sym(f),findsym(sym(f),0.524648)+. subs(sym(f),findsym(sym(f),-0.524648)+. 0.100000*(subs(sym(f),findsym(sym(f),1.650680)+. subs(sym(f),findsym(sym(f),-1.650680); case 5, I=0.945309*subs(sym(f),findsym(sym(f),0)+. 0.393619*(subs(sym(f),findsym(sym(f),0.958572)+. subs(sym(f),findsym(sym(f),-0.958572)+. 0.199532*(subs(sym(f),findsym(sym(f),2.020183)+. subs(sym(f),findsym(sym(f),-2.020183); end11. IntQBXF1求第一类切比雪夫积分function q = IntQBXF1(func,n)format long;pi = 3.1415926535;q = 0;A = zeros(n,1);x = zeros(n,1);for i=1:n A(i) = pi/n; x(i) = cos(pi*(2*i-1)/2/n); y(i) = subs(sym(func), findsym(sym(func), x(i); q = q + A(i)*y(i);end12. IntQBXF2求第二类切比雪夫积分function q = IntQBXF2(func,n)format long;pi = 3.1415926535;q = 0;A = zeros(n,1);x = zeros(n,1);for i=1:n A(i) = sin(i*pi)/(n+1)*sin(i*pi)/(n+1)*pi/(n+1); x(i) = cos(pi*i/(n+1); y(i) = subs(sym(func), findsym(sym(func), x(i); q = q + A(i)*y(i);end13. DblTraprl用梯形公式求重积分function q=DblTraprl(f,a,A,b,B,m,n)if(m=1 & n=1) %梯形公式 q=(B-b)*(A-a)/4)*(subs(sym(f),findsym(sym(f),a,b)+. subs(sym(f),findsym(sym(f),a,B)+. subs(sym(f),findsym(sym(f),A,b)+. subs(sym(f),findsym(sym(f),A,B); else %复合梯形公式 C=4*ones(n+1,m+1); C(1,:)=2; C(:,1)=2; C(n+1,:)=2; C(:,m+1)=2; C(1,1)=1; C(1,m+1)=1; C(n+1,1)=1; C(n+1,m+1)=1; %C矩阵endF=zeros(n+1,m+1);q=0;for i=0:n for j=0:m x=a+i*(A-a)/n; y=b+j*(B-b)/m; F(i+1,j+1)=subs(sym(f),findsym(sym(f),x,y); q=q+F(i+1,j+1)*C(i+1,j+1); endendq=(B-b)*(A-a)/4/m/n)*q;14. DblSimpson用辛普森公式求重积分function q=DblSimpson(f,a,A,b,B,m,n)if(m=1 & n=1) %辛普森公式 q=(B-b)*(A-a)/9)*(subs(sym(f),findsym(sym(f),a,b)+. subs(sym(f),findsym(sym(f),a,B)+. subs(sym(f),findsym(sym(f),A,b)+. subs(sym(f),findsym(sym(f),A,B)+. 4*subs(sym(f),findsym(sym(f),(A-a)/2,b)+. 4*subs(sym(f),findsym(sym(f),(A-a)/2,B)+. 4*subs(sym(f),findsym(sym(f),a,(B-b)/2)+. 4*subs(sym(f),findsym(sym(f),A,(B-b)/2)+. 16*subs(sym(f),findsym(sym(f),(A-a)/2,(B-b)/2); else %复合辛普森公式 q=0; for i=0:n-1 for j=0:m-1 x=a+2*i*(A-a)/2/n; y=b+2*j*(B-b)/2/m; x1=a+(2*i+1)*(A-