Home > voicebox > hypergeom1f1.m

hypergeom1f1

PURPOSE ^

HYPERGEOM1F1 Confluent hypergeometric function, 1F1 a.k.a Kummer's M function [h,l]=(a,b,z,tol,maxj)

SYNOPSIS ^

function [h,l,alg]=hypergeom1f1(a,b,z,tol,maxj,th)

DESCRIPTION ^

 HYPERGEOM1F1 Confluent hypergeometric function, 1F1 a.k.a Kummer's M function [h,l]=(a,b,z,tol,maxj)

  Inputs: a,b,z  input arguments
                 a and b must be real scalars, z can be a real matrix
                 The function is undefined if b is a non-positive integer
          tol    Optional tolerance [default 1e-10]
          maxj   Optional iteration limit [default 500]
          th    Threshold for choice of algorithm [default 30]

 Outputs: h      output result (size of z): 1F1(a;b;z) = M(a;b;z)
          l      actual number of iterations (size of z)
          alg     algorithm used

 This routine is functionally equivalent to the symbolic toolbox function hypergeom(a,b,z) but much faster.
 The function calculated is the solution to z M'' + (b-z) M' - a M = 0 where
 M' and M'' are 1st and 2nd derivatives of M with respect to z.
 This routine is closely based on taylorb1f1.m from [4] which is explained in Sec 3.2 of [3].

 Special cases and relations [2]:
        M(a;b;z) = Inf for integer b<=0
        M(a;b;0) = 1
        M(0;b;z) = 1
        M(a;a;z) = exp(z)
        M(1;2;2z) = exp(z)*sinh(z)/z
        M(a;b;z) = exp(z)*M(b-a;b;-z)    (13.2.39 from [2])
        M(a;a+1;-z) = exp(z)*M(1;a+1;z) = a*gamma(z,a)/z^a
        M(a;b;z) = M(a-1;b;z) + M(a;b+1;z)*z/b
        (b-a)M(a-1,b,z)+(2a-b+z)M(a,b,z)-aM(a+1,b,z)=0     [1] 13.4.1
        b(b-1)M(a,b-1,z)+b(1-b-z)M(a,b,z)+z(b-a)(a,b+1,z)=0     [1] 13.4.2

 References:
 [1]    M. Abramowitz and I. A. Stegun, editors.
       Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.
       Dover, New York, 9 edition, 1964. ISBN 0-486-61272-4.
 [2]    F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.
       NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8.
 [3]    J. Pearson. Computation of hypergeometric functions. Masterís thesis, Oxford University, 2009.
       URL http://¨people.maths.ox.ac.uk/¨porterm/¨research/-pearson_final.pdf.
 [4]    J. Pearson. Computation of hypergeometric functions. Matlab code, Oxford University, 2009.
       URL http://¨people.maths.ox.ac.uk/¨porterm/¨research/-hypergeometricpackage.zip.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [h,l,alg]=hypergeom1f1(a,b,z,tol,maxj,th)
0002 % HYPERGEOM1F1 Confluent hypergeometric function, 1F1 a.k.a Kummer's M function [h,l]=(a,b,z,tol,maxj)
0003 %
0004 %  Inputs: a,b,z  input arguments
0005 %                 a and b must be real scalars, z can be a real matrix
0006 %                 The function is undefined if b is a non-positive integer
0007 %          tol    Optional tolerance [default 1e-10]
0008 %          maxj   Optional iteration limit [default 500]
0009 %          th    Threshold for choice of algorithm [default 30]
0010 %
0011 % Outputs: h      output result (size of z): 1F1(a;b;z) = M(a;b;z)
0012 %          l      actual number of iterations (size of z)
0013 %          alg     algorithm used
0014 %
0015 % This routine is functionally equivalent to the symbolic toolbox function hypergeom(a,b,z) but much faster.
0016 % The function calculated is the solution to z M'' + (b-z) M' - a M = 0 where
0017 % M' and M'' are 1st and 2nd derivatives of M with respect to z.
0018 % This routine is closely based on taylorb1f1.m from [4] which is explained in Sec 3.2 of [3].
0019 %
0020 % Special cases and relations [2]:
0021 %        M(a;b;z) = Inf for integer b<=0
0022 %        M(a;b;0) = 1
0023 %        M(0;b;z) = 1
0024 %        M(a;a;z) = exp(z)
0025 %        M(1;2;2z) = exp(z)*sinh(z)/z
0026 %        M(a;b;z) = exp(z)*M(b-a;b;-z)    (13.2.39 from [2])
0027 %        M(a;a+1;-z) = exp(z)*M(1;a+1;z) = a*gamma(z,a)/z^a
0028 %        M(a;b;z) = M(a-1;b;z) + M(a;b+1;z)*z/b
0029 %        (b-a)M(a-1,b,z)+(2a-b+z)M(a,b,z)-aM(a+1,b,z)=0     [1] 13.4.1
0030 %        b(b-1)M(a,b-1,z)+b(1-b-z)M(a,b,z)+z(b-a)(a,b+1,z)=0     [1] 13.4.2
0031 %
0032 % References:
0033 % [1]    M. Abramowitz and I. A. Stegun, editors.
0034 %       Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.
0035 %       Dover, New York, 9 edition, 1964. ISBN 0-486-61272-4.
0036 % [2]    F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.
0037 %       NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8.
0038 % [3]    J. Pearson. Computation of hypergeometric functions. Masterís thesis, Oxford University, 2009.
0039 %       URL http://¨people.maths.ox.ac.uk/¨porterm/¨research/-pearson_final.pdf.
0040 % [4]    J. Pearson. Computation of hypergeometric functions. Matlab code, Oxford University, 2009.
0041 %       URL http://¨people.maths.ox.ac.uk/¨porterm/¨research/-hypergeometricpackage.zip.
0042 %
0043 
0044 %      Copyright (C) Mike Brookes 2016
0045 %      Version: $Id: hypergeom1f1.m 9980 2017-07-09 17:43:08Z dmb $
0046 %
0047 %   VOICEBOX is a MATLAB toolbox for speech processing.
0048 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0049 %
0050 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0051 %   This program is free software; you can redistribute it and/or modify
0052 %   it under the terms of the GNU General Public License as published by
0053 %   the Free Software Foundation; either version 2 of the License, or
0054 %   (at your option) any later version.
0055 %
0056 %   This program is distributed in the hope that it will be useful,
0057 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0058 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0059 %   GNU General Public License for more details.
0060 %
0061 %   You can obtain a copy of the GNU General Public License from
0062 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0063 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0065 if nargin<4 || isempty(tol)
0066     tol=1e-10;
0067 end
0068 if nargin<5 || isempty(maxj)
0069     maxj=500;
0070 end
0071 if nargin<6 || isempty(th)
0072     th=30;
0073 end
0074 a1=a-1; % define some useful constants
0075 b1=b-1;
0076 ba=b-a;
0077 h=zeros(size(z));
0078 l=zeros(size(z));
0079 nz=numel(z);
0080 for i=1:nz
0081     y=z(i);
0082     q=0;  % break criterion initially false
0083     if abs(y)<th            % small |y|
0084         % compute using the series from 13.2.2 in [2]
0085         % sum(j=0:inf,prod(a+(0:j-1)) y^j /(prod(b+(0:j-1)) j!))
0086         alg=1;
0087         d=1;
0088         g=1;
0089         jlim=0;             % find j beyond which terms definitely decrease
0090         k=(b1+y)^2-4*(a1)*y;
0091         if k>=0
0092             jlim=max(jlim,0.5*(sqrt(k)-(b1+y)));
0093         end
0094         k=(b1-y)^2+4*(a1)*y;
0095         if k>=0
0096             jlim=max(jlim,0.5*(sqrt(k)-(b1-y)));
0097         end
0098         jlim=min(maxj,jlim);
0099         for j=1:maxj
0100             d=d*y*(a1+j)/(j*(b1+j));  % d = A(j)-A(j-1) = (A(j-1)-A(j-2))*r(j)*z from [4]
0101             g=g+d;                      % g = A(j) from [4]
0102             p=abs(d)<tol*abs(g);
0103             if q && p && j>=jlim
0104                 break
0105             end
0106             q=p;
0107         end
0108     elseif y>0 % big positive y
0109         % see 13.7.1 combined with 13.2.3 [2]. Valid for large y unless a is a non-positive integer
0110         alg=2;
0111         d=1;
0112         g=1;
0113         jlim=1;             % find j beyond which terms definitely increase
0114         k=(ba-1-a+y)^2+4*a*(ba-1);
0115         if k>=0
0116             jlim=max(jlim,0.5*(sqrt(k)-(ba-a-1+y)));
0117         end
0118         k=(ba-a-1-y)^2+4*a*(ba-1);
0119         if k>=0
0120             jlim=max(jlim,0.5*(sqrt(k)-(ba-a-1-y)));
0121         end
0122         jlim=min(maxj,jlim);
0123         for j=1:jlim
0124             d=d*(ba-1+j)*(j-a)/(j*y);
0125             g=g+d;
0126             p=abs(d)<tol*abs(g);
0127             if q && p
0128                 break
0129             end
0130             q=p;
0131         end
0132         [gl,gs]=gammalns([a b]);
0133         g=gs(1)*gs(2)*exp(y+gl(2)-gl(1)+(a-b)*log(y))*g;
0134     else   % big negative y
0135         % Combine M(a;b;z) = exp(z)*M(b-a;b;-z) (13.2.39 from [2]) with algorithm 2
0136         alg=3;
0137         d=1;
0138         g=1;
0139                 jlim=1;             % find j beyond which terms definitely increase
0140         k=(a1-ba+y)^2+4*a1*ba;
0141         if k>=0
0142             jlim=max(jlim,0.5*(sqrt(k)-(a1-ba+y)));
0143         end
0144         k=(a1-ba-y)^2+4*a1*ba;
0145         if k>=0
0146             jlim=max(jlim,0.5*(sqrt(k)-(a1-ba-y)));
0147         end
0148         jlim=min(maxj,jlim);
0149         for j=1:jlim
0150             d=d*(a1+j)*(ba-j)/(j*y);
0151             g=g+d;
0152             p=abs(d)<tol*abs(g);
0153             if q && p
0154                 break
0155             end
0156             q=p;
0157         end
0158         [gl,gs]=gammalns([ba b]);
0159         g=gs(1)*gs(2)*exp(gl(2)-gl(1)-a*log(-y))*g;
0160     end
0161     h(i)=g;
0162     l(i)=j;
0163 end

Generated on Tue 19-Sep-2017 12:07:31 by m2html © 2003