Matlab实现——Adaptive Quadrature Using Simpson’s Rule

欢迎前往个人博客 驽马点滴 和视频空间 哔哩哔哩-《挨踢日志》

%Input - f is the integrand input as a string ¡¯f¡¯
%    - a and b are upper and lower limits of integration
%   - tol is the tolerance
% Output  - SRmat is the table of values
%                - quad is the quadrature value
%                - err is the error estimate
%Initialize values
SRmat = zeros(30, 6);
iterating=0;
done=1;
SRvec=zeros(1, 6);
SRvec=srule(f, a, b, tol);
SRmat(1, 1:6)=SRvec;
m=1;
state=iterating;
while(state==iterating)
     n=m;
     for j=n:-1:1
         p=j;
         SR0vec=SRmat(p,:);
         err=SR0vec(5);
         tol=SR0vec(6);
if (tol<=err)
      %Bisect interval,apply Simpson¡¯s rule
      %recursively, and determine error
       state=done;
       SR1vec=SR0vec; SR2vec=SR0vec;
       a=SR0vec(1); b=SR0vec(2); c=(a+b)/2;
       err=SR0vec(5); tol=SR0vec(6); tol2=tol/2;
       SR1vec=srule(f,a,c,tol2);
       SR2vec=srule(f,c,b,tol2);
       err=abs(SR0vec(3)-SR1vec(3)-SR2vec(3))/10;
%Accuracy test
if (err<tol)
     SRmat(p, :)=SR0vec;
     SRmat(p,4)=SR1vec(3)+SR2vec(3);  SRmat(p,5)=err;
else
     SRmat(p+1:m+1,:)=SRmat(p:m,:);
     m=m+1;
     SRmat(p,:)=SR1vec;  SRmat(p+1,:)=SR2vec;
     state=iterating;
end
end
     end
end
quad=sum(SRmat(:,4));
err=sum(abs(SRmat(:,5)));
SRmat=SRmat(1:m,1:6);

 

srule.m

 

%Program7.5(Simpson's Rule)

function Z=srule(f,a,b,tol)

%Input- f is the integrand input as a string ??f??

%- a and b are upper and lower limits of integration

%- tol is the tolerance

% Output - Z is a 1x6 vector [a b < xmlnamespace prefix ="st1" ns ="urn:schemas-microsoft-com:office:smarttags" />S S2 err tol1]

h=(b-a)/2;

C=zeros(1,3);

C=feval(f,[a (a+b)/2 b]);

S=h*(C(1)+4*C(2)+C(3))/3;

S2=S;

tol1=tol;

err=tol;

Z=[a b S S2 err tol1];

 

 

 

fun.m

 

function f=fun(x)

f=sin(4*x).*exp(-2*x);

 

 

 

Untitled.m

 

a=0;

b=3;

tol=0.00001;

f='fun';

[SRmat, quad, err]=adapt(f, a, b, tol)

欢迎前往个人博客  驽马点滴 和视频空间 哔哩哔哩-《挨踢日志》