Matlab数值积分程序集合.doc

Matlab数值积分程序集合[图书馆+网络收集] 近来学习数值积分,手头积累了不少程序,也拿来和各位朋友分享一下。。。主要是来自数值积分教材和网络，基本的原理也就不打算多说了，随便搜索一下就可以得到，那就开始上代码了，呵呵，非原创，但是全部验证过，有疑问可以给我e-mail： 1 梯形数值积分的MATLAB主程序 function T=rctrap(fun,a,b,m)??????? %fun 函数，a 积分上限 b积分下限 m 递归次数 n=1;h=b-a; T=zeros(1,m+1); x=a; T(1)=h*(feval(fun,a)+feval(fun,b))/2; for i=1:m ?????????? h=h/2; n=2*n; s=0; ????????? for k=1:n/2 ?????????? x=a+h*(2*k-1); s=s+feval(fun,x); end T(i+1)=T(i)/2+h*s; end T=T(1:m);e.g 运行后屏幕显示 精确值Fs，用rctrap计算的递归值T和T与精确值Fs的绝对误差wT ) exp((-x^.2./2)./(sqrt(2*pi)))T=rctrap(fun,0,pi/2,14), syms tfi=int(exp((-t^2)/2)/(sqrt(2*pi)),t,0, pi/2);Fs= double(fi), wT= double(abs(fi-T)) fun = ??? @(x)exp((-x^.2./2)./(sqrt(2*pi))) T = Columns 1 through 7 ??? 1.4168??? 1.3578??? 1.3313??? 1.3195??? 1.3142??? 1.3119??? 1.3109 Columns 8 through 14 ??? 1.3105??? 1.3103??? 1.3102??? 1.3102??? 1.3101??? 1.3101??? 1.3101 Fs = ??? 0.4419 wT = Columns 1 through 7 ??? 0.9749??? 0.9159??? 0.8894??? 0.8776??? 0.8723??? 0.8700??? 0.8690 Columns 8 through 14 ??? 0.8686??? 0.8684??? 0.8683??? 0.8683??? 0.8683??? 0.8682??? 0.8682 2 复合辛普森（Simpson）数值积分的MATLAB主程序 function y=comsimpson(fun,a,b,n) % fun 函数 a 积分上限 b积分下限 n 分割小区间数 z1=feval (fun,a)+ feval (fun,b);m=n/2; 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 (fun,x2); end for k=3:2:2*m x3=x+k*h; z3= z3+4*feval (fun,x3); end y=(z1+z2+z3)*h/3; 由于Matlab自带了 quad就是这个算法 所以比较少自己编 3 龙贝格数值积分的MATLAB主程序 function [RT,R,wugu,h]=romberg(fun,a,b, wucha,m)%fun被积函数 a，b积分上下限 wucha两次相邻迭代绝对差值 m 龙贝格积分表最大行数%RT 龙贝格积分表 R 数值积分结果 wucha 误差估计 h 最小步长 n=1;h=b-a; wugu=1; x=a;k=0; RT=zeros(4,4); RT(1,1)=h*(feval(fun,a)+feval(fun,b))/2; while((wuguwucha)(km)|(k4)) ????????? k=k+1;?? h=h/2; s=0; ?????????? for j=1:n ???????????? x=a+h*(2*j-1); s=s+feval(fun,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 wugu=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,wugu,h]=romberg(F,0,1.5,1.e-8,13)syms