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)

## CROSS-REFERENCE INFORMATION

This function calls:
 error Display message and abort function. ipermute Inverse permute array dimensions. num2str Convert number to string. (Fast version) permute Permute array dimensions. shiftdim Shift dimensions. squeeze Remove singleton dimensions.
This function is called by:
 Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics dir2enc Transform a dir. spectrum to an encountered. (Used in spec2spec) dir2k1d Transforms directional spectrum to wavenumber spectrum. dspec2char Evaluates directional spectral characteristics dspec2dcov Return covariance function given a directional spectrum fourier Returns Fourier coefficients. getcrossspectra Compute the cross spectra by integration imlm Iterated maximum likelihood method for estimating the directional distribution jonswap Calculates (and plots) a JONSWAP spectral density memfun private function for MEM seasim Spectral simulation of a Gaussian sea, 2D (x,t) or 3D (x,y,t) spec2bw Evaluates some spectral bandwidth and irregularity factors spec2char Evaluates spectral characteristics and their covariance spec2mom Calculates spectral moments from spectrum spec2spec Transforms between different types of spectra torsethaugen Calculates a double peaked (swell + wind) spectrum wspecplot Plot a spectral density

## 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 %
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