HYPGF Hypergeometric function F(a,b,c,x) CALL: [y ,abserr] = hypgf(a,b,c,x,abseps,releps) y = F(a,b,c,x) abserr = absolute error estimate a,b,c,x = input parameters abseps = requested absolute error releps = requested relative error HYPGF calculates one solution to Gauss's hypergeometric differential equation: x*(1-x)Y''(x)+[c-(a+b+1)*x]*Y'(x)-a*b*Y(x) = 0 where F(a,b,c,x) = Y1(x) = 1 + a*b*x/c + a*(a+1)*b*(b+1)*x^2/(c*(c+1))+.... Many elementary functions are special cases of F(a,b,c,x): 1/(1-x) = F(1,1,1,x) = F(1,b,b,x) = F(a,1,a,x) (1+x)^n = F(-n,b,b,-x) atan(x) = x*F(.5,1,1.5,-x^2) asin(x) = x*F(.5,.5,1.5,x^2) log(x) = x*F(1,1,2,-x) log(1+x)-log(1-x) = 2*x*F(.5,1,1.5,x^2) NOTE: only real x, abs(x) < 1 and c~=0,-1,-2,... are allowed. Examples: x = linspace(-.99,.99)'; [Sn1,err1] = hypgf(1,1,1,x); plot(x,abs(Sn1-1./(1-x)),'b',x,err1,'r'),set(gca,'yscale','log') [Sn2,err2] = hypgf(.5,.5,1.5,x.^2); plot(x,abs(x.*Sn2-asin(x)),'b',x,abs(x.*err2),'r'),set(gca,'yscale','log')
Check if all input arguments are either scalar or of common size. | |
Create cell array. | |
Display message and abort function. | |
Gamma function. | |
Logarithm of gamma function. | |
Not-a-Number. | |
Write formatted data to string. | |
Compare first N characters of strings ignoring case. | |
Display warning message; disable or enable warning messages. |
Parameter estimates for 2D Weibull data. | |
Mean and variance for the 2D Weibull distribution |
001 function [y ,abserr] = hypgf(a,b,c,x,varargin) 002 %HYPGF Hypergeometric function F(a,b,c,x) 003 % 004 % CALL: [y ,abserr] = hypgf(a,b,c,x,abseps,releps) 005 % 006 % y = F(a,b,c,x) 007 % abserr = absolute error estimate 008 % a,b,c,x = input parameters 009 % abseps = requested absolute error 010 % releps = requested relative error 011 % 012 % HYPGF calculates one solution to Gauss's hypergeometric differential 013 % equation: 014 % 015 % x*(1-x)Y''(x)+[c-(a+b+1)*x]*Y'(x)-a*b*Y(x) = 0 016 % where 017 % F(a,b,c,x) = Y1(x) = 1 + a*b*x/c + a*(a+1)*b*(b+1)*x^2/(c*(c+1))+.... 018 % 019 % 020 % Many elementary functions are special cases of F(a,b,c,x): 021 % 1/(1-x) = F(1,1,1,x) = F(1,b,b,x) = F(a,1,a,x) 022 % (1+x)^n = F(-n,b,b,-x) 023 % atan(x) = x*F(.5,1,1.5,-x^2) 024 % asin(x) = x*F(.5,.5,1.5,x^2) 025 % log(x) = x*F(1,1,2,-x) 026 % log(1+x)-log(1-x) = 2*x*F(.5,1,1.5,x^2) 027 % 028 % NOTE: only real x, abs(x) < 1 and c~=0,-1,-2,... are allowed. 029 % 030 % Examples: 031 % x = linspace(-.99,.99)'; 032 % [Sn1,err1] = hypgf(1,1,1,x); 033 % plot(x,abs(Sn1-1./(1-x)),'b',x,err1,'r'),set(gca,'yscale','log') 034 % [Sn2,err2] = hypgf(.5,.5,1.5,x.^2); 035 % plot(x,abs(x.*Sn2-asin(x)),'b',x,abs(x.*err2),'r'),set(gca,'yscale','log') 036 % 037 038 039 040 % Reference: 041 % Kreyszig, Erwin (1988) 042 % Advanced engineering mathematics 043 % John Wiley & Sons, sixth edition, pp 204. 044 045 % Revised pab 30.08.2004 046 % -added dea3 for more stable error estimation 047 % revised pab 18.04.2001 048 % - updated help header 049 % - added example 050 % - added reference 051 % By Per A. Brodtkorb 17.11.98 052 [errorcode, x, a, b, c] = comnsize(x,a,b ,c); 053 if errorcode > 0 054 error('Requires non-scalar arguments to match in size.'); 055 end 056 057 [y,abserr] = gethgf(a,b,c,x,varargin{:}); 058 059 return 060 061 function [fsum, err]=gethgf(a,b,c,x,absEps,relEps,Kmax) 062 error(nargchk(4,7,nargin)) 063 if (nargin<7|isempty(Kmax)) 064 Kmax = 10000; 065 end 066 Kmin = 2; 067 if (nargin<6 | isempty(relEps)), 068 relEps = eps*10; 069 end 070 if (nargin<5|isempty(absEps)) 071 absEps = 0;%eps*10; 072 end 073 useDEA = 1; 074 %abseps = max(abseps,eps*10); 075 076 fsum = zeros(size(x)); 077 delta = fsum; 078 err = fsum; 079 080 ok = ~((round(c)==c & c<=0) | abs(x)>1); 081 k1=find(~ok); 082 if any(k1), 083 warning('Illegal input: c = 0,-1,-2,... or abs(x)>1') 084 fsum(k1) = NaN; 085 err(k1) = NaN; 086 end 087 k0=find(ok & abs(x)==1); 088 if any(k0) 089 cmab = c(k0)-a(k0)-b(k0); 090 fsum(k0) = exp(gammaln(c(k0))+gammaln(cmab)-... 091 gammaln(c(k0)-a(k0))-gammaln(c(k0)-b(k0))); 092 err(k0) = eps; 093 k00 = find(real(cmab)<=0); 094 if any(k00) 095 err(k0(k00)) = nan; 096 fsum(k0(k00)) = nan; 097 end 098 end 099 k=find(ok & abs(x)<1); 100 if any(k), 101 delta(k) = ones(size(k)); 102 fsum(k) = delta(k); 103 104 k1 = k; 105 E = cell(1,3); 106 E{3} = fsum(k); 107 converge = 'n'; 108 for ix=0:Kmax-1, 109 delta(k1) = delta(k1).*((a(k1)+ix)./(ix+1)).*((b(k1)+ix)./(c(k1)+ ... 110 ix)).*x(k1); 111 fsum(k1) = fsum(k1)+delta(k1); 112 113 E(1:2) = E(2:3); 114 E{3} = fsum(k1); 115 116 if ix>Kmin 117 if useDEA, 118 [Sn, err(k1)] = dea3(E{:}); 119 k00 = find((abs(err(k1))) <= max(absEps,abs(relEps.*fsum(k1)))); 120 if any(k00) 121 fsum(k1(k00)) = Sn(k00); 122 end 123 if (ix==Kmax-1) 124 fsum(k1) = Sn; 125 end 126 k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.*fsum(k1))))); 127 if any(k0),% compute more terms 128 nk=length(k0);%# of values we have to compute again 129 E{2} = E{2}(k0); 130 E{3} = E{3}(k0); 131 else 132 converge='y'; 133 break; 134 end 135 else 136 err(k1) = 10*abs(delta(k1)); 137 k0 = (find((abs(err(k1))) > max(absEps,abs(relEps.* ... 138 fsum(k1))))); 139 if any(k0),% compute more terms 140 nk=length(k0);%# of values we have to compute again 141 else 142 converge='y'; 143 break; 144 end 145 end 146 k1 = k1(k0); 147 end 148 end 149 if ~strncmpi(converge,'y',1) 150 disp(sprintf('#%d values did not converge',length(k1))) 151 end 152 end 153 %ix 154 return 155 function [result,abserr] =dea3(E0,E1,E2) 156 %DEA3 Extrapolate a slowly convergent sequence 157 % 158 % CALL: [result,abserr] =dea3(E0,E1,E2) 159 % 160 % result = extrapolated value 161 % abserr = absolute error estimate 162 % E0,E1,E2 = 3 values of a convergent sequence to extrapolate 163 % 164 % DEA3 attempts to extrapolate nonlinearly to a better estimate of the 165 % sequence's limiting value, thus improving the rate of 166 % convergence. Routine is based on the epsilon algorithm 167 % of P. Wynn. 168 % 169 % Example % integrate sin(x) from 0 to pi/2 170 % for I = 5:7 171 % NPARTS = 2.^I; 172 % x = linspace(0,pi/2,2.^I+1); 173 % Ei(I-4) = trapz(x,sin(x)) 174 % end 175 % [En, err] = dea3(Ei(1),Ei(2),Ei(3)) 176 % truErr = Ei-1 177 178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% -*- Mode: Matlab -*- %%%%%%%%%%%%%%%%%%%%%%%%%%%% 179 %% dea3.m --- 180 %% Author : Per Andreas Brodtkorb 181 %% Created On : Wed Aug 18 12:40:23 2004 182 %% Last Modified By: Per Andreas Brodtkorb 183 %% Last Modified On: Wed Aug 18 12:51:49 2004 184 %% Update Count : 3 185 %% Status : Unknown, Use with caution! 186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 187 188 % REFERENCES "Acceleration de la convergence en analyse numerique", 189 % C. Brezinski, "Lecture Notes in Math.", vol. 584, 190 % Springer-Verlag, New York, 1977. 191 192 193 error(nargchk(3,3,nargin)) 194 ten = 10.0d0; 195 one = 1.0d0; 196 small = eps; %spacing(one) 197 delta2 = E2 - E1; 198 delta1 = E1 - E0; 199 err2 = abs(delta2); 200 err1 = abs(delta1); 201 tol2 = max(abs(E2),abs(E1)) * small; 202 tol1 = max(abs(E1),abs(E0)) * small; 203 204 result = zeros(size(E0)); 205 abserr = result; 206 207 k0 = find( ( err1 <= tol1 ) | (err2 <= tol2)); 208 if any(k0) 209 %C IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE 210 %C ACCURACY, CONVERGENCE IS ASSUMED. 211 result(k0) = E2(k0); 212 abserr(k0) = err1(k0) + err2(k0) + E2(k0)*small*ten; 213 end 214 k1 = find(~(( err1 <= tol1 ) | (err2 <= tol2))); 215 if any(k1) 216 ss = one./delta2(k1) - one./delta1(k1); 217 k2 = k1(find((abs(ss.*E1(k1)) <= 1.0d-3))); 218 if any(k2) 219 result(k2) = E2(k2); 220 abserr(k2) = err1(k2) + err2(k2) + E2(k2)*small*ten; 221 end 222 k3 = find(~(abs(ss.*E1(k1)) <= 1.0d-3)); 223 if any(k3) 224 k33 = k1(k3); 225 result(k33) = E1(k33) + one./ss(k3); 226 abserr(k33) = err1(k33) + err2(k33) + abs(result(k33)-E2(k33)); 227 end 228 end 229 return 230 % end subroutine dea3 231 232 233 function [fsum, delta]=gethgf2(a,b,c,z) 234 Kmax=10000; fac=1; tmps=0; 235 tmp = (gamma(a).*gamma(b))./(fac.*gamma(c)); 236 delta=tmp; 237 tol=eps; 238 for ix=1:Kmax 239 fac = fac*ix; 240 tmp = (gamma(a+ix).*gamma(b+ix))./(fac.*gamma(c+ix)); 241 delta = tmp*(z.^ix); 242 243 % don't bother with more terms than we use 244 if (max(abs(delta)) <= tol) , break ,end 245 % update 246 tmps = tmps + delta 247 248 end 249 if ix>=Kmax, disp('Some values did not converge'), end 250 fsum=(gamma(c)./(gamma(a).*gamma(b))).*tmps; 251 delta=(gamma(c)./(gamma(a).*gamma(b))).*delta; 252 return 253 254 255 256 257 258 259 260
Comments or corrections to the WAFO group