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

para la función Gamma, con

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,

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:

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

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

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

realiza una secuencia de operaciones básicas para obtener

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

Referencias

  1. Pugh thesis
Este artículo ha sido escrito por Wikipedia. El texto está disponible bajo la licencia Creative Commons - Atribución - CompartirIgual. Pueden aplicarse cláusulas adicionales a los archivos multimedia.