Гамма-функция и ее свойства

Content

DSPL-2.0 is free digital signal processing algorithms library

Distributed under v3 LGPL v3 license

GitHub project page.

Found a mistake? Select it with the mouse and press ctrl+enter
Введение

Гамма функция находит очень широкое применение в прикладном анализе. С гамма-функцией связана функции Бесселя используемые при синтезе фильтров и спектральном анализе, а также другие специальные функции: бета-функция, К-функции, G-функции. В статистике широко используется гамма-распределение, частными случаями которого являются экспоненциальное распределение и распределение хи-квадрат. В данной статье введено понятие гамма-функции, приведены ее основные свойства, а также показан алгоритм ее численного расчета.

Определение гамма-функции

В математике вводится понятие факториала для натурального числа:

(1)
При этом можно заметить, что

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

(3)
Для натуральных значений аргумента гамма функция совпадает со значением факториала:

(4)
При этом для любых комплексных значений справедливо равенство:

(5)
Данное рекуррентное соотношение является очень важным и используется при расчете гамма-функции. Приведем также формулу дополнения:

(6)
Можно заметить, что при отрицательных значениях , , при этом гамма-функция для отрицательного аргумента может быть вычислена по формуле:

(7)
Необходимо отметить, что при целых ,

и гамма-функция претерпевает разрыв. График гамма-функции для вещественного аргумента представлен на рисунке 1.




Рисунок 1: График гамма-функции вещественного аргумента

Некоторые значения гамма-функции

Рассмотрим некоторые значения гамма-функции. Из выражения (4) следует, что:

(8)
Рассмотрим , для этого воспользуемся выражением (6):

(9)
Рассмотрим

, для этого воспользуемся выражением (5):

(10)
Рассмотрим

, для этого воспользуемся выражением (7):

(11)

Расчет гамма-функции

Теперь рассмотрим очень важный вопрос, касающийся расчета гамма-функции. Для этого рассмотрим несколько возможных интервалов аргумента .



Пусть , тогда в соответствии с (5) можно вычислить:

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

Пусть , тогда можно снова воспользоваться выражением (5), которое можно переписать к виду:

(13)
При этом , и если продолжать, то можно свести к вычислению гамма-функции в интервале .

Рассмотрим на примере:

(14)
То есть опять свели к вычислению гамма-функции в интервале .

Пусть теперь , тогда при вычислении по формуле (7) можно рекуррентно вычислять путем сведения к гамма-функции в интервале .

Теперь для вычисления гамма-функции необходимо получить алгоритм ее расчета при . На практике для этого производят аппроксимацию гамма функции на данном интервале в виде:


(15)
где и — полиномы 8 степени:

(16)
Коэффициенты полиномов аппроксимации подобраны так, чтобы обеспечивать наименьшую ошибку аппроксимации. Значения коэффициентов полиномов приведены в таблице:
1 2 3 4 5 6 7 8
6.65e+4 -3.61e+4 -3.14e+4 866.97 629.33 -379.8 24.77 -1.716
-1.15e+5 -1.35e+5 4.76e+3 2.25e+4 -3107.8 -1015.2 315.35 -30.84

Таким образом, используя полиномиальную аппроксимацию и рекуррентные соотношения можно вычислить значения гамма-функции для любого вещественного аргумента. Программная реализация функции расчета гамма-функции на C приведена ниже.

При численном расчете гамма-функции необходимо соблюдать осторожность, так как скорость роста гамма-функции как у факториала и при 32 битной разрядности процессора вычислять гамма-функцию без переполнения разрядности можно только для аргумента меньшего 170. Например, значение гамма-функции для аргумента 50 равно 6e+62.

Выводы

Таким образом в данной статье мы ввели понятие гамма-функции, рассмотрели ее свойства и привели алгоритм численного расчета гамма-функции на основе полиномиальной аппроксимации. В конце приведен пример программной реализации гамма-функции

Программная реализациия гамма-функции

Следующая программа использует рекуррентные соотношения для расчета гамма-функции.


#include <stdio.h>
#include <stdlib.h>
#define _USE_MATH_DEFINES
#include <math.h>


//***************************************************
// аппроксимация гамма-функции в интервале от 1 до 2
// отношением полиномов 8 степени
double gammaapprox(double x)
{
  double p[]={-1.71618513886549492533811e+0,
               2.47656508055759199108314e+1,
              -3.79804256470945635097577e+2,
               6.29331155312818442661052e+2,
               8.66966202790413211295064e+2,
              -3.14512729688483675254357e+4,
              -3.61444134186911729807069e+4,
               6.64561438202405440627855e+4};

  double q[]={-3.08402300119738975254353e+1,
                3.15350626979604161529144e+2,
               -1.01515636749021914166146e+3,
               -3.10777167157231109440444e+3,
                2.25381184209801510330112e+4,
                4.75584627752788110767815e+3,
               -1.34659959864969306392456e+5,
               -1.15132259675553483497211e+5};
  double z = x - 1.0;
  double a = 0.0;
  double b = 1.0;
  int i;
  for(i = 0; i < 8; i++)
  { 
    a =(a + p[i]) * z;
    b = b * z + q[i];
  }
  return (a / b + 1.0);
}


//***************************************************
// Гамма-функция вещественного агрумента
// возвращает значение гамма-функции аргумента z
double gamma(double z)
{
  
  // рекурентное соотношение для 0
  if((z>0.0)&&(z<1.0))   
    return gamma(z+1.0)/z;   
  
  // рекурентное соотношение для z>2
  if(z>2)
    return (z-1)*gamma(z-1); 

  // рекурентное соотношение для z<=0
  if(z<=0)
    return M_PI/(sin(M_PI*z)*gamma(1-z));

  // 1<=z<=2 использовать аппроксимацию
  return gammaapprox(z); 
}


//***************************************************
// Основная программа для рассчета значения 
// гамма-функции вещественного аргумента
int main(){
  float z = 12.0;
  double g = gamma((double)z);
  printf("gamma(%.2f) = %e\n", z,g);
  return 0;
}

Page update: 24.07.2020 (14:58:42)