Приближение Ланцоша

В математике приближение Ланцоша — это метод численного вычисления гамма-функции , опубликованный Корнелиусом Ланцошем в 1964 году. Это практическая альтернатива более популярному приближению Стирлинга для вычисления гамма-функции с фиксированной точностью.

Введение

Приближение Ланцоша состоит из формулы

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

для гамма-функции, с

$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 .$

Здесь g — действительная константа , которую можно выбирать произвольно с учетом ограничения Re( z + g + 1/2 0 . ) > Коэффициенты p , зависящие от g , вычислить несколько сложнее (см. ниже). Хотя формула, изложенная здесь, действительна только для аргументов в правой комплексной полуплоскости , ее можно распространить на всю комплексную плоскость с помощью формулы отражения :

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

Ряд A сходится и может быть усечен для получения приближения с желаемой точностью. При выборе подходящего g потребуется всего около 5–10 членов ряда (обычно небольшого целого числа) для вычисления гамма-функции с типичной одинарной или двойной точностью с плавающей запятой . фиксированное g Если выбрано , коэффициенты можно рассчитать заранее и благодаря разложению на частные дроби сумма преобразуется в следующий вид:

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

Таким образом, вычисление гамма-функции становится вопросом вычисления лишь небольшого количества элементарных функций и умножения на сохраненные константы. Приближение Ланцоша было популяризировано Numerical Recipes , согласно которому вычисление гамма-функции становится «не намного сложнее, чем другие встроенные функции, которые мы принимаем как должное, такие как sin x или e». ^х Метод также реализован в GNU Scientific Library , Boost , CPython и musl .

Коэффициенты

Коэффициенты имеют вид

$p_{k}(g)={\frac {\sqrt {2\,}}{\pi }}\sum _{\ell =0}^{k}C_{2k+1,\,2\ell +1}\left(\ell -{\tfrac {1}{2}}\right)!{\left(\ell +g+{\tfrac {1}{2}}\right)}^{-(\ell +1/2)}e^{\ell +g+1/2}$

где $C_{n,m}$ представляет ( n , m )-й элемент матрицы коэффициентов полиномов Чебышева , который можно вычислить рекурсивно из этих тождеств:

${\begin{aligned}C_{1,\,1}&=1\\[5px]C_{2,\,2}&=1\\[5px]C_{n+1,\,1}&=-\,C_{n-1,\,1}&{\text{ for }}n&=2,3,4\,\dots \\[5px]C_{n+1,\,n+1}&=2\,C_{n,\,n}&{\text{ for }}n&=2,3,4\,\dots \\[5px]C_{n+1,\,m+1}&=2\,C_{n,\,m}-C_{n-1,\,m+1}&{\text{ for }}n&>m=1,2,3\,\dots \end{aligned}}$

Годфри (2001) описывает, как получить коэффициенты, а также значение усеченного ряда A как матричного произведения .

Вывод

Ланцош вывел формулу из Леонарда Эйлера . интеграла

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

выполнение последовательности основных манипуляций для получения

$\Gamma (z+1)=(z+g+1)^{z+1}e^{-(z+g+1)}\int _{0}^{e}{\Big (}v(1-\log v){\Big )}^{z}v^{g}\,dv,$

и вывод ряда для интеграла.

Простая реализация

Следующая реализация на языке программирования Python работает для сложных аргументов и обычно дает 13 правильных десятичных знаков. Обратите внимание, что исключение наименьших коэффициентов (например, в погоне за скоростью) дает совершенно неточные результаты; коэффициенты необходимо пересчитывать с нуля для разложения с меньшим количеством членов.

from cmath import sin, sqrt, pi, exp


"""
The coefficients used in the code are for when g = 7 and n = 9
Here are some other samples

g = 5
n = 5
p = [
    1.0000018972739440364,
    76.180082222642137322,
    -86.505092037054859197,
    24.012898581922685900,
    -1.2296028490285820771
]

g = 5
n = 7
p = [
    1.0000000001900148240,
    76.180091729471463483,
    -86.505320329416767652,
    24.014098240830910490,
    -1.2317395724501553875,
    0.0012086509738661785061,
    -5.3952393849531283785e-6
]

g = 8
n = 12
p = [
    0.9999999999999999298,
    1975.3739023578852322,
    -4397.3823927922428918,
    3462.6328459862717019,
    -1156.9851431631167820,
    154.53815050252775060,
    -6.2536716123689161798,
    0.034642762454736807441,
    -7.4776171974442977377e-7,
    6.3041253821852264261e-8,
    -2.7405717035683877489e-8,
    4.0486948817567609101e-9
]
"""

g = 7
n = 9
p = [
    0.99999999999980993,
    676.5203681218851,
    -1259.1392167224028,
    771.32342877765313,
    -176.61502916214059,
    12.507343278686905,
    -0.13857109526572012,
    9.9843695780195716e-6,
    1.5056327351493116e-7
]

EPSILON = 1e-07
def drop_imag(z):
    if abs(z.imag) <= EPSILON:
        z = z.real
    return z

def gamma(z):
    z = complex(z)
    if z.real < 0.5:
        y = pi / (sin(pi * z) * gamma(1 - z))  # Reflection formula
    else:
        z -= 1
        x = p[0]
        for i in range(1, len(p)):
            x += p[i] / (z + i)
        t = z + g + 0.5
        y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x
    return drop_imag(y)
"""
The above use of the reflection (thus the if-else structure) is necessary, even though 
it may look strange, as it allows to extend the approximation to values of z where 
Re(z) < 0.5, where the Lanczos method is not valid.
"""

print(gamma(1))
print(gamma(5))   
print(gamma(0.5))

См. также

• Приближение Стирлинга
• Приближение Спуджа

Ссылки

• Годфри, Пол (2001). «Реализация гамма-функции Ланцоша» .
• Ланцос, Корнелиус (1964). «Точная аппроксимация гамма-функции». Журнал Общества промышленной и прикладной математики, серия B: Численный анализ . 1 : 86–96. Бибкод : 1964SJNA....1...86L . дои : 10.1137/0701008 . ISSN   0887-459X . JSTOR   2949767 .
• Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, Б.П. (2007), «Раздел 6.1. Гамма-функция» , Численные рецепты: искусство научных вычислений (3-е изд.), Нью-Йорк: Cambridge University Press, ISBN  978-0-521-88068-8
• Пью, Глендон (2004). Анализ гамма-приближения Ланцоша (PDF) (кандидатская диссертация).
• Тот, Виктор (2005). «Программируемые калькуляторы: приближение Ланцоша» .
• Вайсштейн, Эрик В. «Приближение Ланцоша» . Математический мир .