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

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

Назначение

Рациональная Чебышевская аппроксимация вещественной функции на конечном интервале.

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

ia80r_c вычисляет лучшую в смысле Чебышева рациональную аппроксимацию с некоторым весом вещественной непрерывной на отрезке  [a, b] функции  F (X) : 

 R*[PHI(X)]  =  ( P1 + P2* PHI(X) + ... + PN+1* PHI(X)N ) /
                       / ( Q1 + Q2* PHI(X) + ... + QM+1* PHI(X)M ) ,

причем  R* [ PHI (X)] обладает тем свойством, что на ней достигается минимум 

     max  | R[PHI(X)] - F(X) |  /  | G(X) | .
    [A, B]

Минимизация проводится на классе рациональных функций  R [ РНI (Х)], имеющих степенью числителя  N, а степенью знаменателя  M.

W.J.Cody, W.Fraser, J.F.Hart, Rational Chebyshev Approximation Using Linear Equations, Numerishe Mathematik, 12 (4), 1968.

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

    int ia80r_c (R_fp f, R_fp phi, R_fp g, real *a, real *b,
            real *eps, integer *n, integer *m, real *p, real *q, real *rab,
            integer *ierr)

   Параметры
f - имя вещественной подпрограммы - функции вычисления функции  F (X);
phi - имя вещественной подпрограммы - функции вычисления функции  PHI (X);
g - имя вещественной подпрограммы - функции вычисления функции  G (X);
a, b - заданные нижняя и верхняя границы отрезка, на котоpом производится аппроксимация (тип: вещественный);
eps - заданная точность вычисления аппроксимации (тип: вещественный);
n, m - заданные степень числителя и знаменателя (тип: целый);
p - вещественный вектоp длины  n + 1, в который помещаются коэффициенты числителя;
q - вещественный вектоp длины  m + 1, в который помещаются коэффициенты знаменателя;
rab - вещественный вектоp длины  (n + m + 8) * (n + m + 2) , используемый в подпрограмме как рабочий; в  rab (1) помещается вычисленное значение наилучшего приближения;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
ierr=65 - когда заданная точность не может быть достигнута;
ierr=66 - когда алгоритм не сходится; возможно, что знаменатель имеет корень на заданном отрезке;
ierr=67 - когда сходимость не может быть достигнута в пределах 20 итераций;
ierr=68 - когда алгоритм не сходится из - за сингулярности матрицы системы уравнений для определения  P и  Q;
ierr=69 - когда получена хорошая аппроксимация, но не обязательно наилучшая в смысле Чебышева. 
   Версии:  нет 
   Вызываемые подпрограммы
utia10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы ia80r_c. 
   Замечания по использованию
 

Если  ierr = 65, то числитель и знаменатель имеют общие множители. B этом случае значения  n и  m следует уменьшить на 1.

Bо всех других случаях, когда  ierr ≠ 0, следует попытаться или уменьшить степени числителя и знаменателя, или изменить пределы отрезка, или уменьшить точность вычислений.

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

Пусть требуется аппроксимировать функцию  √ x  на  [0.5, 2.]  следующим образом: 

    x / 4.  +  ( p1 + p2* x ) / ( q1 + q2* x ) .

int main(void)
{
    /* Local variables */
    extern int ia80r_c(R_fp, R_fp ,R_fp , float *, float *, float *, int *,
                       int *, float *, float *, float *, int *);
    static int ierr;
    static float a, b;
    extern float f_c(), g_c(), phi_c();
    static int m, n;
    static float p[2], q[2], rab[40];
    static float eps;

    a = .5f;
    b = 2.f;
    eps = .01f;
    n = 1;
    m = 1;
    ia80r_c((R_fp)f_c, (R_fp)phi_c, (R_fp)g_c, &a, &b, &eps, &n, &m, p, q,
            rab, &ierr);

    printf("\n %16.7e %16.7e %16.7e \n",a, b, eps);
    printf("\n %5i %5i \n",n, m);
    printf("\n %16.7e %16.7e \n",p[0], p[1]);
    printf("\n %16.7e %16.7e \n",q[0], q[1]);
    printf("\n %16.7e %16.7e \n",rab[0], rab[1]);
    printf("\n %5i \n",ierr);
    return 0;
} /* main */

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

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

    ret_val = (float)sqrt(*x) - *x / 4.f;
    return ret_val;
} /* f_c */

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

    ret_val = *x;
    return ret_val;
} /* phi_c */

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

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

    ret_val = (float)sqrt(*x);
    return ret_val;
} /* g_c */


Результаты:

       rab (1)  =  -0.000319 ;    ierr  =  0 ;

       p (1)  =  0.232604 ,    p (2)  =  1.32410 ;

       q (1)  =  1. ,    q (2)  =  1.074789