Aproximación de Lanczos

En matemáticas, la aproximación de Lanczos es un método para evaluar la función gamma numéricamente, publicado por Cornelius Lanczos en 1964. Se trata de una alternativa práctica a la más popular aproximación de Stirling para el cálculo de la función gamma con precisión fija.

Introducción

La aproximación de Lanczos aplica la fórmula

${\displaystyle \Gamma (z+1)={\sqrt {2\pi }}{\left(z+g+{\frac {1}{2}}\right)}^{z+{\frac {1}{2}}}e^{-\left(z+g+{\frac {1}{2}}\right)}A_{g}(z)}$

para la función Gamma, con

${\displaystyle A_{g}(z)={\frac {1}{2}}p_{0}(g)+p_{1}(g){\frac {z}{z+1}}+p_{2}(g){\frac {z(z-1)}{(z+1)(z+2)}}+\cdots .}$

Aquí g es una constante que puede ser elegida arbitrariamente sujeta a la restricción de que Re (z)> 1/2.[1]​ Los coeficientes p, dependientes de g, son ligeramente más difícil de calcular (véase más adelante). Aunque la fórmula como se indica aquí solo es válida para los argumentos en el semiplano complejo derecho, se puede extender a todo el plano complejo mediante la fórmula de reflexión,

${\displaystyle \Gamma (1-z)\;\Gamma (z)={\pi \over \sin \pi z}.}$

La serie A es convergente, y se puede truncar para obtener la aproximación con la precisión deseada. Al elegir una apropiada g (pequeño entero), solo son necesarios de entre 5 a 10 términos de la serie para calcular la función gamma con la típica de precisión simple o doble de coma flotante. Si se elige un g fijo, los coeficientes se pueden calcular de antemano y la suma se refunden en la siguiente forma:

${\displaystyle A_{g}(z)=c_{0}+\sum _{k=1}^{N}{\frac {c_{k}}{z+k}}}$

Por lo tanto el cálculo de la función gamma se convierte en una cuestión de la evaluación de solo un pequeño número de funciones elementales y la multiplicación por constantes almacenadas. La aproximación de Lanczos fue popularizado por Numerical Recipes, según la cual el cálculo de la función Gamma se convierte en "no mucho más difícil que otras funciones incorporadas que damos por sentado, como el sen x o e^x. El método también se implementa en el Biblioteca Científica GNU.

Coeficientes

Los coeficientes están dados por

${\displaystyle p_{k}(g)=\sum _{a=0}^{k}C(2k+1,2a+1){\frac {\sqrt {2}}{\pi }}\left(a-{\begin{matrix}{\frac {1}{2}}\end{matrix}}\right)!{\left(a+g+{\begin{matrix}{\frac {1}{2}}\end{matrix}}\right)}^{-\left(a+{\frac {1}{2}}\right)}e^{a+g+{\frac {1}{2}}}}$

donde C (i, j)denota la (i, j) -ésimo elemento del polinomio de Chebyshev. Los coeficientes de la matriz puede calcularse recursivamente a partir de las identidades

${\displaystyle C(1,1)=1\,}$
${\displaystyle C(2,2)=1\,}$
${\displaystyle C(i,1)=-C(i-2,1)\,}$	${\displaystyle i=3,4,\dots \,}$
${\displaystyle C(i,j)=2C(i-1,j-1)\,}$	${\displaystyle i=j=3,4,\dots \,}$
${\displaystyle C(i,j)=2C(i-1,j-1)-C(i-2,j)\,}$	${\displaystyle i>j=2,3,\dots .}$

Paul Godfrey describe cómo obtener los coeficientes y también el valor de la serie truncada A como un producto de matrices.

Derivación

Lanczos deriva la fórmula de la integral de Leonhard Euler

${\displaystyle \Gamma (z+1)=\int _{0}^{\infty }t^{z}\,e^{-t}\,dt,}$

realiza una secuencia de operaciones básicas para obtener

${\displaystyle \Gamma (z+1)=(z+g+1)^{z+1}e^{-(z+g+1)}\int _{0}^{e}[v(1-\log v)]^{z-{\frac {1}{2}}}v^{g}\,dv,}$

y derivar una serie para la integral.

Implementación simple

La siguiente implementación en el lenguaje de programación Python sirve para argumentos complejos y por lo general da 15 decimales correctos:

def gamma(z):             # great function from Wiki, but maybe could use memorization?
    epsilon = 0.0000001
    def withinepsilon(x):
        return abs(x - abs(x)) <= epsilon

    from cmath import sin,sqrt,pi,exp

    p = [ 676.5203681218851,   -1259.1392167224028,  771.32342877765313,
         -176.61502916214059,     12.507343278686905, -0.13857109526572012,
            9.9843695780195716e-6, 1.5056327351493116e-7]
    z = complex(z)

    # Reflection formula
    if z.real < 0.5:
        result = pi / (sin(pi*z) * gamma(1-z))
    else:
        z -= 1
        x = 0.99999999999980993

        for (i, pval) in enumerate(p):
            x += pval/(z+i+1)

        t = z + len(p) - 0.5
        result = sqrt(2*pi) * t**(z+0.5) * exp(-t) * x

    if withinepsilon(result.imag):
        return result.real
    return result

Véase también

• Aproximación de Stirling
• Aproximación de Spouge

Referencias

1. Pugh thesis

• Godfrey, Paul (2001). «Lanczos Implementation of the Gamma Function». 
• Lanczos, Cornelius (1964). «A Precision Approximation of the Gamma Function». SIAM Journal on Numerical Analysis series B (Society for Industrial and Applied Mathematics) 1: 86-96. ISSN 0887-459X. JSTOR 2949767. doi:10.2307/2949767. 
• Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. (2007), «Section 6.1. Gamma Function», Numerical Recipes: The Art of Scientific Computing (3rd edición), New York: Cambridge University Press, ISBN 978-0-521-88068-8 .
• Pugh, Glendon (2004). An analysis of the Lanczos Gamma approximation (Tesis de PhD). 
• Toth, Viktor (2005). «Programmable Calculators: The Lanczos Approximation». 
• Weisstein, Eric W. «Lanczos Approximation». En Weisstein, Eric W, ed. MathWorld (en inglés). Wolfram Research. 

• Datos: Q6483801