原文地址:matlab数值积分程序集合[转载]作者:卫尉蒙恬1 梯形数值积分的matlab主程序 f nction t=rctrap(f n,a,b,m) %f n 函数,a 积分上限 b积分下限 m 递归次数 n=1;h=b-a; t=zeros(1,m+1); x=a; t(1)=h*(feval_r(f n,a)+feval_r(f n,b))2; for i=1:m h=h2; n=2*n; s=0; for k=1:n2 x=a+h*(2*k-1); s=s+feval_r(f n,x); end t(i+1)=t(i)2+h*s; end t=t(1:m); e.g 运行后屏幕显示 精确值fs,用rctrap计算的递归值t和t与精确值fs的绝对误差wt  f n=@(x) exp((-x^.2.2).(sqrt(2*pi t=rctrap(f n,0,pi2,14), syms t fi=int(exp((-t^2)2)(sqrt(2*pi)),t,0, pi2); fs= do ble(fi), wt= do ble(abs(fi-t)) f n = @(x)exp((-x^.2.2).(sqrt(2*pi t = col mns 1 thro gh 7 1.4168 1.3578 1.3313 1.3195 1.3142 1.3119 1.3109 col mns 8 thro gh 14 1.3105 1.3103 1.3102 1.3102 1.3101 1.3101 1.3101 fs = 0.4419 wt = col mns 1 thro gh 7 0.9749 0.9159 0.8894 0.8776 0.8723 0.8700 0.8690 col mns 8 thro gh 14 0.8686 0.8684 0.8683 0.8683 0.8683 0.8682 0.8682  2 复合辛普森(simpson)数值积分的matlab主程序 f nction y=comsimpson(f n,a,b,n) % f n 函数 a 积分上限 b积分下限 n 分割小区间数 z1=feval_r(f n,a)+ feval_r(f n,b);m=n2; h=(b-a)(2*m); x=a; z2=0; z3=0; x2=0; x3=0; for k=2:2:2*m x2=x+k*h; z2= z2+2*feval_r(f n,x2); end for k=3:2:2*m x3=x+k*h; z3= z3+4*feval_r(f n,x3); end y=(z1+z2+z3)*h3; 由于matlab自带了 q ad就是这个算法 所以比较少自己编 3 龙贝格数值积分的matlab主程序 f nction [rt,r,w g ,h]=romberg(f n,a,b, w cha,m) %f n被积函数 a,b积分上下限 w cha两次相邻迭代绝对差值 m 龙贝格积分表最大行数 %rt 龙贝格积分表 r 数值积分结果 w cha 误差估计 h 最小步长 n=1;h=b-a; w g =1; x=a;k=0; rt=zeros(4,4); rt(1,1)=h*(feval_r(f n,a)+feval_r(f n,b))2; while((w g w cha)&(km)|(k4)) k=k+1; h=h2; s=0; for j=1:n x=a+h*(2*j-1); s=s+feval_r(f n,x); end rt(k+1,1)= rt(k,1)2+h*s; n=2*n; for i=1:k rt(k+1,i+1)=((4^i)*rt(k+1,i)-rt(k,i))(4^i-1); end w g =abs(rt(k+1,k)-rt(k+1,k+1)); end r=rt(k+1,k+1); e.g.  f=inline( 1.(1+x) ); [rt,r,w g ,h]=romberg(f,0,1.5,1.e-8,13) syms x fi=int(1(1+x),x,0,1.5); fs=do ble(fi), wr=do ble(abs(fi-r)), wr1= wr - w g rt = 1.0500 0 0 0 0 0 0.9536 0.9214 0 0 0 0 0.9260 0.9168 0.9165 0 0 0 0.9187 0.9163 0.9163 0.9163 0 0 0.9169 0.9163 0.9163 0.9163 0.9163 0 0.9164 0.9163 0.9163 0.9163 0.9163 0.9163 r = 0.9163 w g = 2.9436e-011 h = 0.0469 fs = 0.9163 wr = 9.8007e-011 wr1 = 6.8571e-011  4 复合梯形法f nction [i,step] = combinetraprl(f,a,b,eps) %f 被积函数 %a,b 积分上下限 %eps 精度 %i 积分结果 %step 积分的子区间数 if(nargin ==3) eps=1.0e-4; end n=1; h=(b-a)2; i1=0; i2=(s bs(sym(f),findsym(sym(f)),a)+s bs(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+(h2)*(s bs(sym(f),findsym(sym(f)),x)+s bs(sym(f),findsym(sym(f)),x1)); end end i=i2; step=n;5 辛普森法 f nction [i,step] = intsimpson(f,a,b,type,eps) %type = 1 辛普森公式 %type = 2 辛普森38公式 %type = 3 复合辛普森公式 if(type==3 && nargin==4) eps=1.0e-4; %缺省精度为0.0001 end i=0; switch type case 1, i=((b-a)6)*(s bs(sym(f),findsym(sym(f)),a)+... 4*s bs(sym(f),findsym(sym(f)),(a+b)2)+... s bs(sym(f),findsym(sym(f)),b)); step=1; case 2, i=((b-a)8)*(s bs(sym(f),findsym(sym(f)),a)+... 3*s bs(sym(f),findsym(sym(f)),(2*a+b)3)+ ... 3*s bs(sym(f),findsym(sym(f)),(a+2*b)3)+s bs(sym(f),findsym(sym(f)),b)); step=1; case 3, n=2; h=(b-a)2; i1=0; i2=(s bs(sym(f),findsym(sym(f)),a)+s bs(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+(h6)*(s bs(sym(f),findsym(sym(f)),x)+...