GAUSSQ Numerically evaluates a integral using a Gauss quadrature. The Quadrature integrates a (2m-1)th order polynomial exactly and the integral is of the form b Int (p(x)* Fun(x)) dx a CALL: [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],p1,p2,....) [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,p1,p2,....) [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,beta,p1,p2,....) int = evaluated integral tol = absolute tolerance abs(int-intold) Fun = string containing the name of the function or a directly given expression enclosed in parenthesis. The function may depend on parameters alpha,beta, pi. xlow,xhigh = integration limits reltol = relative tolerance default=1e-3 wfun = weight function 1 p(x)=1 a =-1, b = 1 Legendre (default) 2 p(x)=1/sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev of the first kind 3 p(x)=sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev of the second kind 4 p(x)=sqrt((x-a)/(b-x)), a = 0, b = 1 5 p(x)=1/sqrt(b-x), a = 0, b = 1 6 p(x)=sqrt(b-x), a = 0, b = 1 7 p(x)=(x-a)^alpha*(b-x)^beta a =-1, b = 1 Jacobi alpha, beta >-1 (default alpha=beta=0) 8 p(x)=x^alpha*exp(-x) a = 0, b = inf generalized Laguerre 9 p(x)=exp(-x^2) a =-inf, b = inf Hermite 10 p(x)=1 a =-1, b = 1 Legendre (slower than 1) trace = for non-zero TRACE traces the function evaluations with a point plot of the integrand. gn = number of base points and weight points to start the integration with (default=2) p1,p2,...= coefficients to be passed directly to function Fun: G = Fun(x,p1,p2,...). Note that int is the common size of xlow, xhigh and p1,p2,.... Example:% a) integration of x.^2 from 0 to 2 and from 1 to 4 % b) integration of x^2*exp(-x) from zero to infinity gaussq('(x.^2)',[0 1],[2 4]) % a) gaussq('(1)',0,inf,[1e-3 8],[],2) % b) gaussq('(x.^2)',0,inf,[1e-3 8],[],0) % b) See also qrule, gaussq2d
Display message and abort function. | |
Construct a string from a function handle. | |
Convert integer to string (Fast version). | |
True if object is a given class. | |
Linear plot. | |
Write formatted data to string. |
Brodtkorb (2004) joint (Scf,Hd) CDF of laboratory data. | |
Brodtkorb (2004) joint (Scf,Hd) CDF from Japan Sea. | |
Brodtkorb et.al. (2000) joint (Scf,Hd) CDF from North Sea. | |
Joint 2D CDF computed as int F(X1 | |
Joint (Scf,Hd) CDF for linear waves with a JONSWAP spectrum. | |
Joint (Scf,Hd) CDF for non-linear waves with JONSWAP spectrum. | |
Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum. | |
Joint (Vcf,Hd) CDF for non-linear waves with JONSWAP spectrum. | |
Mean and variance for the MDIST2D distribution. | |
Myrhaug and Kjeldsen (1987) joint (Scf,Hd) CDF. | |
Joint (Scf,Hd) CDF for linear waves with Ochi-Hubble spectra. | |
Joint (Scf,Hd) CDF for linear waves in space with Ochi-Hubble spectra. | |
Tayfun (1981) CDF of breaking limited wave heights | |
Tayfun (1981) PDF of breaking limited wave heights | |
Tayfun (1990) CDF of large wave heights | |
Tayfun (1990) PDF of large wave heights | |
Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) CDF for nonlinear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) CDF for linear waves in space with Torsethaugen spectra. | |
Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra. | |
Joint 2D Weibull cumulative distribution function | |
Mean and variance for the 2D Weibull distribution |
001 function [int, tol1,ix]= gaussq(fun,xlow,xhigh,tol,trace,varargin) 002 %GAUSSQ Numerically evaluates a integral using a Gauss quadrature. 003 % The Quadrature integrates a (2m-1)th order polynomial exactly 004 % and the integral is of the form 005 % b 006 % Int (p(x)* Fun(x)) dx 007 % a 008 % 009 % CALL: 010 % [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],p1,p2,....) 011 % [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,p1,p2,....) 012 % [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,beta,p1,p2,....) 013 % 014 % int = evaluated integral 015 % tol = absolute tolerance abs(int-intold) 016 % Fun = string containing the name of the function or a directly given 017 % expression enclosed in parenthesis. The function may depend on 018 % parameters alpha,beta, pi. 019 %xlow,xhigh = integration limits 020 % reltol = relative tolerance default=1e-3 021 % wfun = weight function 022 % 1 p(x)=1 a =-1, b = 1 Legendre (default) 023 % 2 p(x)=1/sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev of the 024 % first kind 025 % 3 p(x)=sqrt((x-a)*(b-x)), a =-1, b = 1 Chebyshev of the 026 % second kind 027 % 4 p(x)=sqrt((x-a)/(b-x)), a = 0, b = 1 028 % 5 p(x)=1/sqrt(b-x), a = 0, b = 1 029 % 6 p(x)=sqrt(b-x), a = 0, b = 1 030 % 7 p(x)=(x-a)^alpha*(b-x)^beta a =-1, b = 1 Jacobi 031 % alpha, beta >-1 (default alpha=beta=0) 032 % 8 p(x)=x^alpha*exp(-x) a = 0, b = inf generalized Laguerre 033 % 9 p(x)=exp(-x^2) a =-inf, b = inf Hermite 034 % 10 p(x)=1 a =-1, b = 1 Legendre (slower than 1) 035 % 036 % trace = for non-zero TRACE traces the function evaluations 037 % with a point plot of the integrand. 038 % gn = number of base points and weight points to start the 039 % integration with (default=2) 040 %p1,p2,...= coefficients to be passed directly to function Fun: 041 % G = Fun(x,p1,p2,...). 042 % 043 % Note that int is the common size of xlow, xhigh and p1,p2,.... 044 % Example:% a) integration of x.^2 from 0 to 2 and from 1 to 4 045 % % b) integration of x^2*exp(-x) from zero to infinity 046 % 047 % gaussq('(x.^2)',[0 1],[2 4]) % a) 048 % gaussq('(1)',0,inf,[1e-3 8],[],2) % b) 049 % gaussq('(x.^2)',0,inf,[1e-3 8],[],0) % b) 050 % 051 % See also qrule, gaussq2d 052 053 054 % References 055 % 056 % [1] Golub, G. H. and Welsch, J. H. (1969) 057 % 'Calculation of Gaussian Quadrature Rules' 058 % Mathematics of Computation, vol 23,page 221-230, 059 % 060 % [2] Davis and Rabinowitz (1975) 'Methods of Numerical Integration', page 365, 061 % Academic Press. 062 % 063 % [3]. Stroud and Secrest (1966), 'gaussian quadrature formulas', 064 % prentice-hall, Englewood cliffs, n.j. 065 % 066 % [4] Abromowitz and Stegun (1954) 'Handbook of mathematical functions' 067 068 069 070 071 % tested on: Matlab 5.3 072 % history: 073 % Revised pab 22Nov2004 074 % -Added the possibility of using a function handle. 075 % Revised pab 09.09.2002 076 % -added the possibility of using a inline function. 077 % revised pab 27.03.2000 078 % - fixed a bug for p1,p2,... changed to varargin in input 079 % revised pab 19.09.1999 080 % documentation 081 % by Per A. Brodtkorb 30.03.99, 17.02.99 : 082 % wfun 1: from grule.m in NIT toolbox, see ref [2] 083 % wfun 2-6: see ref [4] 084 % wfun 7-10: Adapted from Netlib routine gaussq.f see ref [1,3] 085 % -accept multiple integrationlimits, int is the common size of xlow and xhigh 086 % -optimized by only computing the integrals which did not converge. 087 % -enabled the integration of directly given functions enclosed in 088 % parenthesis. Example: integration from 0 to 2 and from 2 to 4 for x is done by: 089 % gaussq('(x.^2)',[0 2],[2 4]) 090 091 global ALPHA1 BETA1 ALPHA2 092 wfun=1; 093 if nargin<4| isempty(tol), 094 tol=1e-3; 095 elseif length(tol)==2, 096 wfun=tol(2); 097 tol=tol(1); 098 end 099 100 P0=varargin; 101 NP=length(P0); 102 103 istart = 0; 104 alpha1 = 0; 105 beta1 = 0; 106 FindWeights = 1; 107 108 x = []; 109 y = []; 110 111 switch wfun 112 case 7, 113 istart=2; 114 if ((NP>=1) & (~isempty(P0{1}))), alpha1 = P0{1}; end 115 if ((NP>=2) & (~isempty(P0{1}))), beta1 = P0{2}; end 116 if isempty(ALPHA1)|isempty(BETA1), 117 elseif ALPHA1==alpha1 & BETA1==beta1, 118 FindWeights=0; 119 end 120 ALPHA1=alpha1;BETA1=beta1; 121 %remember what type of weights are saved as global constants 122 case 8, 123 istart=1; 124 if ((NP>=1) & (~isempty(P0{1}))), alpha1 = P0{1}; end 125 if isempty(ALPHA2), 126 elseif ALPHA2==alpha1, 127 FindWeights=0; 128 end 129 ALPHA2=alpha1; 130 %remember what type of weights are saved as global constants 131 otherwise,FindWeights=0; 132 end 133 134 P0(1:istart)=[]; 135 136 gn=2; 137 if nargin <5|isempty(trace) , 138 trace = 0; 139 elseif length(trace)==2, 140 gn = trace(2); 141 trace = trace(1); 142 end 143 144 145 if prod(size(xlow))==1,% make sure the integration limits have correct size 146 xlow=xlow(ones(size(xhigh)));; 147 elseif prod(size(xhigh))==1, 148 xhigh=xhigh(ones(size(xlow)));; 149 elseif any( size(xhigh)~=size(xlow) ) 150 error('The input must have equal size!') 151 end 152 [N M]=size(xlow);%remember the size of input 153 154 num_parameters=NP-istart; 155 if num_parameters>0, 156 ptxt = sprintf('p%d,',1:num_parameters); 157 ptxt(end)=[]; % remove ',' 158 eval(sprintf('[%s]=deal(P0{:});',ptxt)); 159 end 160 %Old call 161 %for ix=1:num_parameters, 162 % eval(['p' int2str(ix) '=P0{ix};']); % variable # i 163 %end 164 165 if (isa(fun,'char') & any(fun=='(')), % & any(fun=='x'), 166 exec_string=['y=',fun ';']; %the call function is already setup 167 else 168 %setup string to call the function 169 if isa(fun,'function_handle') 170 fun = func2str(fun); 171 end 172 exec_string=['y=feval(fun,x']; 173 for ix=1:num_parameters, 174 xvar=['p' int2str(ix)]; % variable # i 175 if eval(['isnumeric(' ,xvar,') & length(',xvar,'(:)) ~=1' ]) , 176 if N*M==1, 177 eval(['[N M]=size(', xvar, ');']); 178 elseif eval(['N*M~=prod(size(', xvar, '))']), 179 error('The input must have equal size!') 180 end 181 eval([xvar, '=' ,xvar ,'(:);']); %make sure it is a column 182 exec_string=[exec_string,',' xvar '(k,ones(1,gn) )']; %enable integration with multiple arguments 183 else 184 exec_string=[exec_string,',' xvar]; 185 end 186 end 187 exec_string=[exec_string,');']; 188 end 189 190 191 nk = N*M;% # of integrals we have to compute 192 k = (1:nk)'; 193 int = zeros(nk,1); 194 tol1 = int; 195 196 wtxt = sprintf('%d_%d',gn,wfun); % number of weights and 197 % weightfunction used 198 cbtxt = sprintf('cb%s',wtxt); %base points string 199 cwtxt = sprintf('cw%s',wtxt); %weights string 200 eval(sprintf('global %s %s ;',cbtxt,cwtxt)); 201 if isempty(eval(['cb',wtxt]))|FindWeights , 202 % calculate the weights if necessary 203 eval(sprintf('[%s,%s]=qrule(gn,wfun,alpha1,beta1);',cbtxt,cwtxt)); 204 end 205 206 %setup mapping parameters and execution strings 207 xlow=xlow(:); 208 jacob=(xhigh(:)-xlow(:))/2; 209 210 switch wfun,% this is clumsy and can written more tidy 211 case {1 ,10}, 212 calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 213 int_string=['(ones(nk,1),:).*y,2).*jacob(k);']; 214 case 2, 215 calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 216 int_string=['(ones(nk,1),:).*y,2);']; 217 case 3, 218 calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 219 int_string=['(ones(nk,1),:).*y,2).*jacob(k).^2;']; 220 case 4, 221 calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));']; 222 int_string=['(ones(nk,1),:).*y,2).*jacob(k)*2;']; 223 case 5, 224 calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));']; 225 int_string=['(ones(nk,1),:).*y,2).*sqrt(jacob(k)*2);']; 226 case 6, 227 calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));']; 228 int_string=['(ones(nk,1),:).*y,2).*sqrt(jacob(k)*2).^3;']; 229 case 7, 230 calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 231 int_string=['(ones(nk,1),:).*y,2).*jacob(k).^(alpha1+beta1+1);']; 232 case {8,9} 233 calcx_string=['(ones(nk,1),:));']; 234 int_string=['(ones(nk,1),:).*y,2);']; 235 otherwise error('unknown option') 236 end 237 238 239 eval(['x=(',cbtxt, calcx_string]); % calculate the x values 240 eval(exec_string); % calculate function values y=fun(x) 241 eval(['int(k)=sum(',cwtxt, int_string]); % do the integration 242 int_old=int; 243 244 if trace==1, 245 x_trace=x(:); 246 y_trace=y(:); 247 end 248 249 % Break out of the iteration loop for three reasons: 250 % 1) the last update is very small (compared to int and compared to tol) 251 % 2) There are more than 11 iterations. This should NEVER happen. 252 253 converge='n'; 254 for i=1:10, 255 gn=gn*2;% double the # of weights 256 wtxt = sprintf('%d_%d',gn,wfun); % # of weights and weight function used 257 %eval(['global cb',wtxt,' cw',wtxt]); 258 cbtxt = sprintf('cb%s',wtxt); %base points string 259 cwtxt = sprintf('cw%s',wtxt); %weights string 260 eval(sprintf('global %s %s ;',cbtxt,cwtxt)); 261 if isempty(eval(cbtxt))|FindWeights , 262 % calculate the weights if necessary 263 eval(sprintf('[%s,%s]=qrule(gn,wfun,alpha1,beta1);',cbtxt,cwtxt)); 264 end 265 266 eval(['x=(',cbtxt, calcx_string]); % calculate the x values 267 eval(exec_string); % calculate function values y=fun(x) 268 eval(['int(k)=sum(',cwtxt, int_string]);% do the integration 269 270 if trace==1, 271 x_trace=[x_trace;x(:)]; 272 y_trace=[y_trace;y(:)]; 273 end 274 275 tol1(k) = abs(int_old(k)-int(k)); %absolute tolerance 276 k = find(tol1 > abs(tol*int)); %| tol1 > abs(tol));%indices to integrals which did not converge 277 278 if any(k),% compute integrals again 279 nk=length(k);%# of integrals we have to compute again 280 else 281 converge='y'; 282 break; 283 end 284 int_old=int; 285 end 286 287 int=reshape(int,N,M); % make sure int is the same size as the integration limits 288 if nargout>1, 289 tol1=reshape(tol1,N,M); 290 end 291 292 if converge=='n' 293 if nk>1 294 if (nk==N*M), 295 disp('All integrals did not converge--singularities likely') 296 else 297 disp(sprintf('%d integrals did not converge--singularities likely', ... 298 nk)) 299 end 300 else 301 disp('Integral did not converge--singularity likely') 302 end 303 end 304 305 if trace==1, 306 plot(x_trace,y_trace,'+') 307 end 308 309 310
Comments or corrections to the WAFO group