【数值分析实验MATLAB】数值积分：梯形公式、辛普森公式、复化梯形公式、复化辛普森公式、龙贝格算法、自适应求积方法

全栈程序员-站长 • 2026年3月19日上午11:32 • 未分类 • 阅读 2

【数值分析实验MATLAB】数值积分：梯形公式、辛普森公式、复化梯形公式、复化辛普森公式、龙贝格算法、自适应求积方法
（1）Matlab代码:
f1:

if x==0 y=1; else %y=exp(-x^2); y=(sin(x))/x; end

梯形公式：

a=0; b=1; tx=(b-a)/2*(f1(a)+f1(b)); tx=vpa(tx,5); tx

simpson公式：

a=0; b=1; xps=(b-a)/6*(f1(a)+4*f1((a+b)/2)+f1(b)); xps=vpa(xps,5); xps

a=0; b=1; e=1e-6; m=1; h=b-a; T(1)=h/2*(f1(a)+f1(b)); while 1 h=h/2; sum=0; for j=1:2^(m-1) sum=sum+f1(a+(2*j-1)*h); end S(1)=T(1)/2+h*sum; for j=1:m S(j+1)=S(j)+(S(j)-T(j))/(4^j-1); end if(abs(S(m+1)-T(m))<=e) fprintf('求得最终解为：%f\n',S(m+1)); fprintf('迭代次数为:%d\n',m); return end T=S; m=m+1; end

复化辛普森公式：

a=0; b=1; e=1e-6; N=20; n=0; h0=(b-a)/2; S0=h0/3*(f1(a)+f1(b)+4*f1((a+b)/2)); while n<N n=n+1; h1=h0/2; sum1=0; sum2=0; for k=1:2^(n-1) sum1=sum1+f1(a+(4*k-2)*h1); end for k=1:2^n sum2=sum2+f1(a+(2*k-1)*h1); end S1=S0/2+h1/3*(-2*sum1+4*sum2); if(abs(S1-S0)<e) fprintf('求得最终解为：%f\n',S1); fprintf('迭代次数为:%d\n',n); return end S0=S1; h0=h1; end

龙贝格算法：

a=0; b=1; e=1e-6; h=b-a; T(1)=h.*(f1(a)+f1(b))./2; m=1; while(1) h=h/2; sum=0; for j=1:2^(m-1) sum=sum+f1(a+h*(2*j-1)); end S(1)=T(1)/2+h*sum; for j=1:m S(j+1)=S(j)+(S(j)-T(j))/(4^j-1); end if(abs(S(m+1)-T(m))<=e) break; end T=S; m=m+1; end fprintf('求得解为：%f\n',S(m+1)); fprintf('迭代次数：%d\n',m);

自适应求积方法：

format short a=0; b=1; e=1e-6; p=[a,b]; p0=p; ep=[e]; m=0; q=0; I=0; while(1) n1=length(ep); n=length(p0); if n==1 break; end h=p0(2)-p0(1); s0=h/6*(f1(p0(1))+4*f1(p0(1)+h/2)+f1(p0(1)+h)); s1=h/12*(f1(p0(1))+4*f1(p0(1)+h/4)+f1(p0(1)+h/2)); s2=h/12*(f1(p0(1)+h/2)+4*f1(p0(1)+h*3/4)+f1(p0(1)+h)); if abs(s0-s1-s2)<=15*ep(1) I=I+s1+s2; p0=p0(2:n); if n1>=2; ep=ep(2:n1); end q=q+1; else m=m+1; p0=[p0(1),p0(1)+h/2,p0(2:n)]; if n1==1 ep=[ep(1)/2,ep(1)/2]; else ep=[ep(1)/2,ep(1)/2,ep(2:n1)]; end if q==0 p=p0; else p=[p(1:q),p0]; end end end fprintf('求得的解为：%f\n',I); fprintf('分层数为%d\n',m); p=p'; disp('分点为');p