Home > wafo > misc > hypgf.m

# hypgf

## PURPOSE

Hypergeometric function F(a,b,c,x)

## SYNOPSIS

[y ,abserr] = hypgf(a,b,c,x,varargin)

## DESCRIPTION

 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')

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. cell Create cell array. error Display message and abort function. gamma Gamma function. gammaln Logarithm of gamma function. nan Not-a-Number. sprintf Write formatted data to string. strncmpi Compare first N characters of strings ignoring case. warning Display warning message; disable or enable warning messages.
This function is called by:
 weib2dfit Parameter estimates for 2D Weibull data. weib2dstat Mean and variance for the 2D Weibull distribution

## SOURCE CODE

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

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