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: This function is called by:

SUBFUNCTIONS ^

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