simpson.m

很多matlab的源代码

function [int,t0,n] = simpson(F,l,tol,kmax)
% SIMPSON Area under string function of t using Simpson's rule.
%
%	[I,ATOL,AN] = SIMPSON(F,T,TOL,N) Area under F(t), a string function
%	An example of a string function is : F(t) = '2*tri(t-1)'
%	T = [Tlo,Thi] = time limits.
%	NOTE: If T is a scalar, the routine uses T = [0 T].
%	TOL = required tolerence [Default: TOL = 1e-6]
%	N =  maximum number of iterations [Default: 12].
%	I = area,  ATOL = actual tolerance in ATOL, AN = evaluation count.
%
%	SIMPSON (with no input agruments) invokes the following example:
%
%	% How much energy is contained in the central lobe of x(t) = sinc(t)
%	% Since E = area(x^2),  we find the area under sinc2(t) from -1 to 1
%         >>[e,atol] = simpson('sinc2(t)',[-1 1])


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998



if nargin==0,help simpson,disp('Strike a key to see results of the example')
pause,[e,atol]=simpson('sinc2(t)',[-1 1]),return,end

if nargin<4;kmax=12;end,
if nargin<3,tol=1e-6;end,
int0=0;n=16;k=0;t0=1;
if length(l)==1,a=0;b=l;else,a=l(1);b=l(2);end,
N=1;h=(b-a)/n;t=a+h*(0:2:n);
y0=eval(F);
s0=sum(y0)+sum(y0(2:n/2));
while t0>tol & k<=kmax,k=k+1;

if k<11,
t=a+h*(1:2:n-1);y=eval(F);s1=sum(y);else,s1=0;N=2*N;n0=n/N;
for j=1:N,
t=a+h*((j-1)*n0+1:2:j*n0-1);
y=eval(F);s1=s1+sum(y);
end,
end

s=s0+2*s1;int=h*(s+2*s1)/3;h=h/2;
s0=s;t0=abs((int-int0)/int);int0=int;n=2*n;

end,

n=n/2;
if k >kmax & t0>tol,disp('WARNING: Recursion limit reached'),end