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
Display message and abort function. | |
Inverse permute array dimensions. | |
Convert number to string. (Fast version) | |
Permute array dimensions. | |
Shift dimensions. | |
Remove singleton dimensions. |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Transform a dir. spectrum to an encountered. (Used in spec2spec) | |
Transforms directional spectrum to wavenumber spectrum. | |
Evaluates directional spectral characteristics | |
Return covariance function given a directional spectrum | |
Returns Fourier coefficients. | |
Compute the cross spectra by integration | |
Iterated maximum likelihood method for estimating the directional distribution | |
Calculates (and plots) a JONSWAP spectral density | |
private function for MEM | |
Spectral simulation of a Gaussian sea, 2D (x,t) or 3D (x,y,t) | |
Evaluates some spectral bandwidth and irregularity factors | |
Evaluates spectral characteristics and their covariance | |
Calculates spectral moments from spectrum | |
Transforms between different types of spectra | |
Calculates a double peaked (swell + wind) spectrum | |
Plot a spectral density |
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
Comments or corrections to the WAFO group