Matlab数值积分程序集合.docx
- 文档编号:17880605
- 上传时间:2023-08-04
- 格式:DOCX
- 页数:13
- 大小:43.27KB
Matlab数值积分程序集合.docx
《Matlab数值积分程序集合.docx》由会员分享,可在线阅读,更多相关《Matlab数值积分程序集合.docx(13页珍藏版)》请在冰点文库上搜索。
Matlab数值积分程序集合
Matlab数值积分程序集合[图书馆+网络收集]
近来学习数值积分,手头积累了不少程序,也拿来和各位朋友分享一下。
。
。
主要是来自数值积分教材和网络,基本的原理也就不打算多说了,随便搜索一下就可以得到,那就开始上代码了,呵呵,非原创,但是全部验证过,有疑问可以给我e-mail:
1梯形数值积分的MATLAB主程序
functionT=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;
fori=1:
m
h=h/2;n=2*n;s=0;
fork=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),symst
fi=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=
Columns1through7
1.4168 1.3578 1.3313 1.3195 1.3142 1.3119 1.3109
Columns8through14
1.3105 1.3103 1.3102 1.3102 1.3101 1.3101 1.3101
Fs=
0.4419
wT=
Columns1through7
0.9749 0.9159 0.8894 0.8776 0.8723 0.8700 0.8690
Columns8through14
0.8686 0.8684 0.8683 0.8683 0.8683 0.8682 0.8682
>>
2复合辛普森(Simpson)数值积分的MATLAB主程序
functiony=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;
fork=2:
2:
2*m
x2=x+k*h;z2=z2+2*feval(fun,x2);
end
fork=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((wugu>wucha)&(k k=k+1; h=h/2;s=0; forj=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; fori=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) symsx fi=int(1/(1+x),x,0,1.5);Fs=double(fi), wR=double(abs(fi-R)),wR1=wR-wugu 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 wugu= 2.9436e-011 h= 0.0469 Fs= 0.9163 wR= 9.8007e-011 wR1= 6.8571e-011 >> 6复合梯形法 function[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=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; whileabs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; fori=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)); end end I=I2; step=n; 7辛普森法 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.0001 end I=0; switchtype case1, 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; case2, 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; case3, n=2; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; whileabs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; fori=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 8牛顿-科茨法 functionI=NewtonCotes(f,a,b,type) %type=1科茨公式 %type=2牛顿-科茨六点公式 %type=3牛顿-科茨七点公式 I=0; switchtype case1, 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)); case2, 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)); case3, 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)); end 9高斯公式 functionI=IntGauss(f,a,b,n,AK,XK) if(n<5&&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)); end ta=(b-a)/2; tb=(a+b)/2; switchn case0, I=2*ta*subs(sym(f),findsym(sym(f)),tb); case1, I=ta*(subs(sym(f),findsym(sym(f)),ta*0.5773503+tb)+... subs(sym(f),findsym(sym(f)),-ta*0.5773503+tb)); case2, 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)); case3, 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)); case4, 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 10区间逐次分半梯形法 functionq=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矩阵 end F=zeros(n+1,m+1); q=0; fori=0: n forj=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); end end q=((B-b)*(A-a)/4/m/n)*q; 11区间逐次分半布尔法 function[I,step]=DDBuer(f,a,b,eps) if(nargin==3) eps=1.0e-4; end; n=1; h=b-a; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h/2; tol=1; whiletol>eps n=n+1; h=h/4; I1=I2; I2=0; fori=0: (4^(n-1)-1) x=a+h*4*i; x1=x+h; x2=x1+h; x3=x2+h; x4=x3+h; I2=I2+(2*h/45)*(7*subs(sym(f),findsym(sym(f)),x)+... 32*subs(sym(f),findsym(sym(f)),x1)+... 12*subs(sym(f),findsym(sym(f)),x2)+... 32*subs(sym(f),findsym(sym(f)),x3)+... 7*subs(sym(f),findsym(sym(f)),x4)); end tol=abs(I2-I1); end I=I2; step=n; 10自适应求积法 functionI=SmartSimpson(f,a,b,eps) if(nargin==3) eps=1.0e-4; end; e=5*eps; I=SubSmartSimpson(f,a,b,e); functionq=SubSmartSimpson(f,a,b,eps) QA=IntSimpson(f,a,b,1,eps); QLeft=IntSimpson(f,a,(a+b)/2,1,eps); QRight=IntSimpson(f,(a+b)/2,b,1,eps); if(abs(QLeft+QRight-QA)<=eps) q=QA; else q=SubSmartSimpson(f,a,(a+b)/2,eps)+SubSmartSimpson(f,(a+b)/2,b,eps); end e.g >>SmartSimpson('x*sin(x)',0,1) ans= 0.3011 >>
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 数值 积分 程序 集合