Home > wafo > misc > simpson.m

simpson

PURPOSE ^

Numerical integration with the Simpson method

SYNOPSIS ^

[area,epsi,a,b] = simpson(x,f,dim)

DESCRIPTION ^

  SIMPSON Numerical integration with the Simpson method
 
  CALL: [area,epsi,a,b] = simpson(x,f,dim)
        [area,epsi,a,b] = simpson(f,dim)
 
     area = result, approximative integral
     epsi = estimate of the absolute error
     a,b  = integration limits
       x  = argument of function (vector or matrix, see below)
       f  = function (vector or matrix, see below)
      dim = dimension to integrate 
 
  Evaluates an approximation of the integral of the
  vector function F(x) wrt X from A to B using the Simpson method. 
 
  X and F must be vectors of the same length,
  OR X must be a vector and F an array whose first non-singleton
  dimension is length(X), then SIMPSON operates along this dimension.
  OR, if X is an array with the same dimension as F the integration
  is performed column-wise. If DIM is given SIMPSON integrates across
  dimension DIM of F. The length of X must be the same as size(F,DIM)).
 
  Notice that if N = length(F) is an even number,
  SIMPSON is used on the N-1 first values and the Trapezoidal rule 
  on the two last values. If X is not equidistant
  spaced use TRAPZ instead.
 
  Example:%
           x = linspace(0,4,201);
           f = exp(-x.^2);
           [area,epsi] = simpson(x,f)
 
  See also  trapz

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [area,epsi,a,b] = simpson(x,f,dim)
002 % SIMPSON Numerical integration with the Simpson method
003 %
004 % CALL: [area,epsi,a,b] = simpson(x,f,dim)
005 %       [area,epsi,a,b] = simpson(f,dim)
006 %
007 %    area = result, approximative integral
008 %    epsi = estimate of the absolute error
009 %    a,b  = integration limits
010 %      x  = argument of function (vector or matrix, see below)
011 %      f  = function (vector or matrix, see below)
012 %     dim = dimension to integrate 
013 %
014 % Evaluates an approximation of the integral of the
015 % vector function F(x) wrt X from A to B using the Simpson method. 
016 %
017 % X and F must be vectors of the same length,
018 % OR X must be a vector and F an array whose first non-singleton
019 % dimension is length(X), then SIMPSON operates along this dimension.
020 % OR, if X is an array with the same dimension as F the integration
021 % is performed column-wise. If DIM is given SIMPSON integrates across
022 % dimension DIM of F. The length of X must be the same as size(F,DIM)).
023 %
024 % Notice that if N = length(F) is an even number,
025 % SIMPSON is used on the N-1 first values and the Trapezoidal rule 
026 % on the two last values. If X is not equidistant
027 % spaced use TRAPZ instead.
028 %
029 % Example:%
030 %          x = linspace(0,4,201);
031 %          f = exp(-x.^2);
032 %          [area,epsi] = simpson(x,f)
033 %
034 % See also  trapz
035 
036 % Tested on: Matlab 5.3
037 % History:
038 % revised by pab 14.05.2000
039 %  - added dim
040 %  - added trapezoidal rule when lengt(f) is even
041 %  - changed the documentation according to the new changes.
042 % revised by es 28.09.1999 (documentation, help)
043 % by pab 1997, updated 16.06.1998
044 
045 if 0 % Old call
046   if nargin<2,
047     f = x;
048     [f,nshift] = shiftdim(f)
049     [nx m]=size(f);
050     fsiz(1)=1;fsiz(2)=1;
051     if nx<m
052       tmp=m;m=nx;nx=tmp; 
053       f=f.';
054       % fsiz(1)=m;  fsiz(2)=1;
055     else
056       %fsiz(1)=1;fsiz(2)=m;
057     end
058     switch m
059       case 1,  a=1;b=nx; bm1=nx-1;
060       case 2,  a=f(1,1);b=f(nx,1);bm1=f(nx-1,1); f=f(:,2);
061       otherwise, error('Wrong dimension of input! dim must be 2xN, 1xN, Nx2 or Nx1 ') 
062     end
063     m=1;
064   end  
065   
066   if nargin==2,
067     fsiz=size(f);
068     [n m]=size(f);
069     [nx mx]=size(x);
070     
071     %make sure x is a column vector
072     if nx*mx==n | nx*mx==m  , x=x(:);nx=max(nx,mx);end 
073     
074     a=x(1,:);b=x(nx,:);bm1=x(nx-1,:); 
075     if nx~=n
076       tmp=m;m=n;n=tmp; 
077       f=f.';
078       fsiz(2)=1;
079       if (nx~=n) | (ndims(f)>2)
080     error('Fx must have dimensions which equals the  size of x ')
081       end
082     else
083       fsiz(1)=1;
084     end
085   end  
086   
087 
088   if nx<3, 
089     error('The vector must have more than 3 elements!')
090   end
091   
092   if (mod(nx, 2) == 0), % checking if n is even
093     f(nx,:)=[];
094     nx=nx-1; 
095     
096     disp('Warning: Changed the integration limits ')
097     disp(['         from ' num2str(b) ' to ' num2str(bm1) ])
098     disp('         so that the length of f(x) is odd.')
099     
100     b=bm1;
101   end
102   
103   hn=2*(b-a)/(nx-1);
104   
105   
106   % Midpoint approximation
107   Un=hn.*sum(f(2:2:(nx-1),1:m)); 
108   
109   % Trapeziodal approximation
110   Tn=hn/2.*(f(1,1:m) + 2*sum(f(3:2:(nx-2),1:m)) + f(nx,1:m)); 
111   
112   % Simpson's approximation
113   area=(Tn+2*Un)/3; 
114   
115   % Asymptotically the simpson approximation is a better estimate 
116   % than both Tn and Un. If this is the case, 
117   % then an estimate of the absolute error is
118   
119   epsi = abs(Tn-Un)./4; 
120   %fsiz
121   %area
122   area=squeeze(reshape(area,fsiz));
123   epsi=squeeze(reshape(epsi,fsiz));
124   
125   
126 
127 else % NEW CALL
128  
129   perm = []; nshifts = 0;
130   if nargin == 3,                 % simpson(x,f,dim)
131     perm = [dim:max(ndims(f),dim) 1:dim-1];  
132     f = permute(f,perm);
133   elseif nargin==2 & length(f)==1 % simpson(f,dim)
134     dim = f; f = x; x=[]
135     perm = [dim:max(ndims(f),dim) 1:dim-1];
136     f = permute(f,perm);
137   else                            % simpson(x,f) or simpson(f)
138     if nargin < 2, f = x; x=[]; end
139     [f,nshifts] = shiftdim(f);
140   end
141   
142   fsiz = size(f);
143   m = fsiz(1);
144   if m> 5 , n=3;else, n=1;end
145   if isempty(x) % assume unit spacing
146     a=1;b=m;bmn=m-n;
147   else
148     if isempty(perm),
149       x = shiftdim(x);
150       xsiz=size(x);
151     else
152       xsiz = size(x);
153       if prod(xsiz)==prod(fsiz)
154     x = permute(x,perm);
155       end
156     end
157     if xsiz(1) ~= m
158       error('length(x) must equal length of first non-singleton dim of f.');
159     end
160     a =x(1,:);b=x(m,:);
161     bmn=x(m-n,:); 
162   end
163     
164   if m<3, 
165     error('The vector must have more than 3 elements!')
166   end
167   
168   if (mod(m, 2) == 0), % checking if m is even
169     
170     if n==1 % Use trapezoidal  from bmn to b
171       area = (f(m,:)+f(m-1,:))/2.*(b-bmn); 
172     else    % Use 4-points Newton-Cotes rule
173       area = (f(m,:)+3*(f(m-1,:)+f(m-2,:))+ f(m-3,:))/8.*(b-bmn); 
174     end
175     m = m-n;     
176     hn=2*(bmn-a)/(m-1);
177   else
178     area = 0;
179     hn=2*(b-a)/(m-1);
180   end
181   
182   % Midpoint approximation
183   Un = hn.*sum(f(2:2:(m-1),:)); 
184   
185   % Trapeziodal approximation
186   Tn = hn/2.*(f(1,:) + 2*sum(f(3:2:(m-2),:)) + f(m,:)); 
187   
188   % Simpson's approximation
189   area = area+(Tn+2*Un)/3; 
190   
191   % Asymptotically the simpson approximation is a better estimate 
192   % than both Tn and Un. If this is the case, 
193   % then an estimate of the absolute error is
194   if nargout>1,  epsi = abs(Tn-Un)/4; end
195   
196   fsiz(1) = 1;
197   area = reshape(area,[ones(1,nshifts),fsiz]);
198   if ~isempty(perm), area = ipermute(area,perm); end
199 end
200

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003