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

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

Назначение

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

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

Пусть на интервале  [u, v] задана функция  f (x). Вычисляется дробно - рациональная функция вида 

                  p+1                                     q+1
   R (x) = (   ∑   A k * Tk - 1 ( z ) )  /  (   ∑   B k * Tk - 1 ( z ) )  ≡
                  k = 1                                    k = 1
                    p+1                          q+1
           ≡   (   ∑   C k x k - 1 )  /  (   ∑   D k x k - 1 ) ,
                   k = 1                          k = 1

где  z = ( 2x - u - v ) / ( v - u ) для  x ∈ [u , v],  Tk ( z ) - полином Чебышева 1 - го рода степени  k,  B1 = D1 = 1,  p, q - заданные натуральные числа, удовлетворяющая условиям интерполяции: 

           R ( x j ) = f ( x j )

в узлах  x j = 0.5 * ( u + v - ( v - u ) cos y j ), где  y j = ( j - 0.5 ) * π / ( p + q + 1 ),  j = 1, 2, ..., p + q + 1.

Для найденного  R (x) в узлах сетки  x k = 0.5 * ( u + v + ( v - u ) cos f k ),   f k = ( k - 1 ) * π / ( p + q + 1 ),  k = 1, 2, ..., p + q + 2, вычисляются характеристики точности интерполяции 

          E1 = max k | Δ f k | ,           E2 = max k | Δ k / f k | ,
     EPS1 = E1 / min k | Δ k |  ,   EPS2 = E2 / min k | Δ k |  ,
 где
   f k = f ( x k ) ,    Δ k = f ( x k ) - R ( x k ) .

Curtis A., Osborne M.R. The construction of minimax rational approximations to functions, Computer J., 1966, v.9, N 3, 286 - 293.

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

    int ii20r_c (integer *ip, integer *iq, real *u, real *v,
        R_fp funct, real *a, real *b, real *al, real *be, real *rab,
            real *rab1, real *rab2, real *rab3, real *rab4, integer *ierr)

   Параметры
ip, iq - заданные степени числителя  p и знаменателя  q,  ip ≥ 0,  iq ≥ 0 (тип: целый);
u, v - заданные нижняя и верхняя границы интервала интерполяции (тип: вещественный);
funct - имя вещественной подпрограммы - функции вычисления значений  f (x);
a - вещественный вектоp длины  ip + 1 значений  a k = a[k];
b - вещественный вектоp длины  iq + 1 значений  b k = b[k];
al - вещественный вектоp длины  ip + 1 значений  c k = al[k];
be - вещественный вектоp длины  iq + 1 значений  d k = be[k];
rab - вещественный вектоp длины  iq, используемый в подпрограмме как рабочий;
rab1 - вещественный двумерный массив размера  3 * (2 ( ip + iq + 1) + 1), используемый в подпрограмме как pабочий;
rab2 - вещественный двумерный массив размера  6 * ( ip + iq + 2), используемый как рабочий; после окончания работы подпрограммы в нем запоминаются хаpактеристики интерполяции  E1 = rab2 (1, 1),  E2 = rab2 (2, 1),  Eps1 = rab2 (3, 1),  Eps2 = rab2 (4, 1);
rab3 - вещественный вектоp длины  ip + 1, используемый в подпрограмме как рабочий;
rab4 - вещественный двумерный массив размера  iq * iq, используемый в подпрограмме как pабочий;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
ierr=65 - когда вычисленный знаменатель искомого приближения имеет корень в интервале  [u, v];
ierr=66 - когда матрица для определения коэффициентов приближения вырождена. 
   Версии:  нет 
   Вызываемые подпрограммы
uti i10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы i i20r_c. 
   Замечания по использованию:  нет

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

int main(void)
{
    /* Local variables */
    extern int ii20r_c(int *, int *, float *, float *, R_fp, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, int *);
    static int ierr;
    static float a[3], b[3];
    static float u, v;
    extern float funct_c();
    static float be[3], al[3];
    static int ip, iq;
    static float rab[2], rab1[33] /* was [3][11] */,
                         rab2[36] /* was [6][6] */,
                 rab3[3], rab4[4] /* was [2][2] */;

#define rab2_ref(a_1,a_2) rab2[(a_2)*6 + a_1 - 7]

    ip = 2;
    iq = 2;
    u = -1.f;
    v = 1.f;
    ii20r_c(&ip, &iq, &u, &v, (R_fp)funct_c, a, b, al, be, rab, rab1, rab2,
            rab3, rab4, &ierr);

    printf("\n  %5i %5i %16.7e %16.7e %5i \n", ip, iq, u, v, ierr);
    printf("\n  %16.7e %16.7e %16.7e \n", a[0],a[1],a[2]);
    printf("\n  %16.7e %16.7e %16.7e \n", b[0],b[1],b[2]);
    printf("\n  %16.7e %16.7e %16.7e \n", al[0],al[1],al[2]);
    printf("\n  %16.7e %16.7e %16.7e \n", be[0],be[1],be[2]);
    printf("\n  %16.7e %16.7e %16.7e %16.7e \n",
            rab2_ref(1, 1), rab2_ref(2, 1), rab2_ref(3, 1), rab2_ref(4, 1));
    return 0;
} /* main */


float funct_c(float *x)
{
    /* System generated locals */
    float ret_val;

    /* Builtin functions */
    double exp(double);

    ret_val = (float)exp(*x);
    return ret_val;
} /* funct_c */


Результаты:

       a   =  (1.,   .480197,   .0392021) ;
       b   =  (1.,  -.480197,   .0392021) ;

       al  =  (1.,   .499789,   .0816032) ;
       be  =  (1.,  -.499789,   .0816032) ;

       ierr  =  0

       rab2 (1, 1)  =   .00024 ;     rab2 (2, 1)  =   .000089 ;
       rab2 (3, 1)  =  7.3884 ;      rab2 (4, 1)  =  1.05524