Гамма-функция и ее применение

Содержание

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

Определение гамма-функции
В математике вводится понятие факториала для натурального числа:
(1)
При этом можно заметить, что
(2)
Гамма-функция , распространяет понятие факториала на дробные, отрицательные и даже комплексные значения аргумента . Гамма функция не выражается через элементарные функции, но может быть представлена как интеграл вида:
(3)
Для натуральных значений аргумента гамма функция совпадает со значением факториала:
(4)
При этом для любых комплексных значений справедливо равенство:
(5)
Данное рекуррентное соотношение является очень важным и используется при расчете гамма-функции. Приведем также формулу дополнения:
(6)
Можно заметить, что при отрицательных значениях , , при этом гамма-функция для отрицательного аргумента может быть вычислена по формуле:
(7)
Необходимо отметить, что при целых , и гамма-функция претерпевает разрыв. График гамма-функции для вещественного аргумента представлен на рисунке 1.


Для просмотра SVG графики Вам необходимо обновить браузер

Рисунок 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;
        for(int 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){
       if((z>0.0)&&(z<1.0))   return gamma(z+1.0)/z;     // рекурентное соотношение для 0
       if(z>2)                return (z-1)*gamma(z-1);   // рекурентное соотношение для z>2
       if(z<=0)               return M_PI/(sin(M_PI*z)*gamma(1-z)); // рекурентное соотношение для z<=0
       return gammaapprox(z); // 1<=z<=2 использовать аппроксимацию
}


//*************************************************************************************************************

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


Любые вопросы и пожелания вы можете оставить в гостевой книге, на форуме, или прислать по электронной почте admin@dsplib.ru


Система Orphus
Любое копирование материалов сайта без разрешения автора запрещено.
Разработка и дизайн Бахурин Сергей.