What's the formula used for kummerU?

Hi! Thanks for reading!

I need to write a kummer U function for myself.

For kummerU funtion, it is said in the doc that:

But it's only for Re(a)>0 and Re(z)>0, while I want the formula for a<0. So I try to read its source code:

function y = kummerU(a,b,x)
%KUMMERU    Kummer's confluent hypergeometric U function.
%   U = kummerU(a,b,x) computes the value of Kummer's U function
%     U = 1/gamma(a) * int(exp(-x*t)*t^(a-1)*(1+t)^(b-a-1),t,0,inf).
%
%   Reference:
%     [1] Abramowitz, Milton; Stegun, Irene A., eds., "Chapter 13",
%     Handbook of Mathematical Functions with Formulas, Graphs, and
%     Mathematical Tables, New York: Dover, pp. 504, 1965.

%   Copyright 2014 The MathWorks, Inc.

y = sym.useSymForNumeric(@kummerU, a, b, x);
end

There is nothing for me. So does the function

sym.useSymForNumeric(@kummerU, a, b, x)

There is no str "kummerU" inside it.

So how can I find out how matlab realize kummerU with negative a?

Or it would be better if you know the formula for negative a.

Any suggestion is appreciated!

[Walter Roberson]

https://cradpdf.drdc-rddc.gc.ca/PDFS/unc57/p527324.pdf

However there are special cases for negative integer a, in which case the formula becomes a polynomial

[祥宇 崔]

@Walter Roberson

Thanks! It helps!

Answer 1

[Walter Roberson]

Execute

regexprep(char(evalin(symengine, 'expose(kummerU)')) , '\\n', '\n')

to be shown the formula used.

[祥宇 崔]

Sometimes it's just hard to find the function we want. There are too many in matlab.

Thanks, experienced Roberson!

This is exactly what I want!!!!

[祥宇 崔]

@Walter Roberson

It seems that kummerU uses laguerre polynomial:

if domtype(ra) = DOM_INT && 0 <= ra && iszero(frac(a)) && ra < Pref::autoExpansionLimit() then
   return((-1)^ra*ra!*laguerreL(ra, b - 1, z))
end_if; 

Then, for laguerreL:

if domtype(n) = DOM_INT && 0 <= n then
   return(orthpoly::laguerre(n, a, x))
end_if; 

For the function:

orthpoly::laguerre(n, a, x)

I look through the code roughly:

proc(n, a, x)
       name orthpoly::laguerre;
       local i, expand_coeffs, r0, r1;
     begin
       if args(0) <> 3 then
         error(message(\"symbolic:orthpoly:IncorrectNumberOfArguments\"))
       end_if; 
       if contains({n, x}, undefined) then
         return(undefined)
       end_if; 
       if contains({n, x}, RD_NAN) then
         return(RD_NAN)
       end_if; 
       n := specfunc::makeInteger(n); 
       if n = FAIL || type(n) = DOM_INT && n < 0 || n = -infinity then
         error(message(\"symbolic:specfunc:ExpectingNonNegativeIntegerOrSymbol\"))
       end_if; 
       if type(n) = stdlib::Infinity then
         return(undefined)
       end_if; 
       if ~domtype(n) = DOM_INT then
         return(procname(args()))
       end_if; 
       if ~testtype(a, Type::Arithmetical) then
         error(message(\"symbolic:orthpoly:InvalidSecondArgument\"))
       end_if; 
       if contains({DOM_RAT, DOM_INT, DOM_FLOAT, DOM_COMPLEX}, domtype(a)) && contains({DOM_RAT, DOM_INT, DOM_FLOAT, DOM_COMPLEX}, domtype(x)) then
         expand_coeffs := FALSE
       else 
         expand_coeffs := TRUE
       end_if; 
       if ~testtype(x, Type::Arithmetical) then
         error(message(\"symbolic:orthpoly:InvalidThirdArgument\"))
       end_if; 
       if iszero(n + a) then
         develinfo(1, \"using special formula\"); 
         return(((-1)^n/n!)*x^n)
       end_if; 
       if domtype(x) = DOM_IDENT || type(x) = \"_index\" then
         case n
           of 0 do
             return(1)
           of 1 do
             return(1 + a - x)
           otherwise
             r0 := poly(1, [x]); 
             r1 := poly(1 + a - x, [x]); 
             assert(expand_coeffs); 
             for i from 1 to n - 1 do
               [r1, r0] := [mapcoeffs(poly((2*i + 1 + a - x)/(i + 1), [x])*r1 - poly((i + a)/(i + 1), [x])*r0, expand), r1]
             end_for; 
             return(op(r1, 1))
         end_case
       else 
         case n
           of 0 do
             return(1)
           of 1 do
             return(1 + a - x)
           otherwise
             r0 := 1; 
             r1 := expand(1 + a - x); 
             if expand_coeffs then
               for i from 1 to n - 1 do
                 [r1, r0] := [expand(((2*i + 1 + a - x)*r1 - (i + a)*r0)/(i + 1)), r1]
               end_for
             else 
               for i from 1 to n - 1 do
                 [r1, r0] := [((2*i + 1 + a - x)*r1 - (i + a)*r0)/(i + 1), r1]
               end_for
             end_if; 
             return(r1)
         end_case
       end_if
     end_proc'

I think it's based on polynomial summation.

So in conclusion, I think for my question:

Matlab use polynomials instead of integration to realize kummerU function. It makes sense since it's usually easier to calculate polynomial than numerical integration.