Текст подпрограммы и версий
sf74r_c.zip  		Тексты тестовых примеров
tsf74r_c.zip 

Подпрограмма:  sf74r_c

Назначение

Вычисление модифицированных функций Бесселя первого рода   In (x) для последовательности целых индексов.

Математическое описание

sf74r_c вычисляет значения модифицированных функций Бесселя первого рода  Ip (x) для действительного аргумента  x и последовательности целых индексов  p от  n до 0. Расчеты осуществляются по рекуррентной формуле: 

     Ip-1(x) = (2p/x)*Ip(x) + Ip+1(x) .

Ошибка вычисления  Ip (x) не возрастает в pекуppентном процессе, применяемом при убывании  p.

Вычисления начинаются с задания значений  Ik + 1 = 0,   Ik = 1 и числа  k.

Величина  k выбирается в зависимости от значений  x и  n. Затем последовательно вычисляются значения функций для p = k - 1, k - 2, ..., 0. Для достижения заданной точности  eps рекурсия повторяется с большим  k до тех пор, пока разность между pезультатами двух последних рекурсий не станет меньше  eps.

1.  Абрамовиц М., Стиган И. Справочник по специальным функциям - М. Наука, 1979.
2.  Goldstein M., Thaler R.M. Recurrence techniques for the calculation of Bessel functions. - M TAC, 1959, V. 13, N 66.

Использование 

    int sf74r_c (real *x, integer *n1, real *eps, logical *s,
            real *t, real *t1)

   Параметры
x - заданное значение аргумента  x (тип: вещественный);
n1 - максимальный порядок последовательности функций   Ip (x), увеличенный на единицу (тип: целый);
eps - относительная точность вычисления последней функции последовательности  In (x) (тип: вещественный);
s - логическая переменная, задающая режим работы подпрограммы; при этом:
если  s = TRUE_, то вычисляются функции  Ip(x),
если  s = FALSE_, то вычисляются функции e - x Ip(x);
t - вещественный вектоp длины n1 вычисленных значений функций Бесселя порядков от 0 до  n;   Ij (x) = t (j + 1),   j = 0, ..., n  (n = n1 - 1);
t1 - вещественный вектоp длины n1, используемый как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  Для значений  x, при которых  e - x меньше положительного наименьшего представимого в машине числа, значение s полагается равным FALSE_

Пример использования 

int main(void)
{
    /* Initialized data */
    static float x = 5.f;
    static float eps = 1e-9f;
    static int n1 = 21;

    /* Local variables */
    extern int sf74r_c(float *, int *, float *, logical *, float *, float *);
    static logical s;
    static float t[21], t1[21];
    int i;

    s = TRUE_;
    sf74r_c(&x, &n1, &eps, &s, t, t1);

    for (i = 1; i <= 19; i += 3) {
    printf("\n %16.7e %16.7e %16.7e ",t[i-1], t[i], t[i+1]);
    }
    printf("\n\n ");
    return 0;
} /* main */


Результаты:

       t(5)    =  5.1082347636e+00
       t(11)  =  4.580044192e-03
       t(16)  =  1.0479776754e-06
       t(21)  =  5.0242393580e-11