Текст подпрограммы и версий
ia83r_c.zip , ia83d_c.zip 		Тексты тестовых примеров
tia83r_c.zip , tia83d_c.zip

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

Назначение

Построение алгебраического полинома наилучшего равномерного приближения для таблично заданной функции.

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

Пусть в узлах сетки  x1 < x2 < ...< xM заданы значения табличной функции  y :  y1, y2 ..., yM. Требуется построить алгебраический полином 

                         N
         P*(x)  =   ∑   ak xk
                       k =0

заданной степени N ≤ M - 2, который минимизирует максимальное отклонение 

                 h ( y,P )   =    max   | yk - P(x) |   
                                     1≤k≤M

среди всех полиномов P (x) степени N.

Среди узлов  x1, x2, ..., xM существуют точки альтернанса  x1* < x2* < ... < x*N + 2, в которых 

                 yk - P* (xk*)  =  (- 1)k h* , 

 где  h*  =   max    | yk - P*(xk ) | , 
                 1≤k≤M

h* - максимальное отклонение полинома наилучшего приближения P* (x).

Задача решается на основе итерационного R - алгоритма В.Н.Малоземова, требующего задания начального набора N + 2 точек x1 < x2 < ...< xN + 2, каждая из которых должна совпадать с одним из узлов исходной сетки. Подпрограмма находит точки альтернанса  x1* < x2* < ...< x*N + 2, не обязательно совпадающие с начальным набором точек.

Сб. "Вопросы теории и элементы программного обеспечения минимаксных задач", под ред. Демьянова В.Ф., Малоземова В.Н., Л., Изд - во Ленингр. ун - та, 135 - 138, 1977.

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

    int ia83r_c(integer *m, integer *n, real *x, real *y,
        integer *in, integer *iter, real *a, real *ax, integer *jn,
        real *rab, integer *ierr)

   Параметры
m - заданное число узлов сетки (тип: целый);
n - заданная степень полинома, n ≤ m - 2 (тип: целый);
x - вещественный вектор длины m, в котором задаются значения узлов сетки;
y - вещественный вектор длины m, в котором задаются значения табличной функции;
in - целый вектор длины n + 2 номеров начального набора точек; может изменяться самой программой в процессе вычисления;
iter - число выполненных при решении задачи итераций (тип: целый);
a - вещественный вектор длины n + 1 вычисленных значений коэффициентов полинома наилучшего приближения, A ( i ) = ai,  i = 1, 2, ...,n + 1;
ax - вещественный вектор длины n + 2 вычисленных значений точек альтернанса;
jn - целый вектор длины n + 2 номеров вычисленных значений точек альтернанса;
rab - вещественный двумерный массив размера 2 * (n + 2), используемый как рабочий; после окончания работы подпрограммы rab (1, 1) = h*;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом
ierr=65 - когда m < n + 2;
ierr=66 - когда номера начального набора точек не упорядочены строго по возрастанию; или упорядоченность нарушается в процессе работы программы;
ierr=67 - когда значения узлов сетки  xi,  i = 1, 2, ..., m не упорядочены строго по возрастанию.
ierr=68 - сходимость не может быть достигнута в пределах 30 итераций. 
   Версии
ia83d_c - построение алгебраического полинома наилучшего равномерного приближения для таблично заданной функции в режиме вычислений с удвоенной точностью. 
   Вызываемые подпрограммы
utia10_c - подпрограмма выдачи диагностических сообщений при работе ia83r_c;
utia11_c - подпрограмма выдачи диагностических сообщений при работе ia83d_c; 
   Замечания по использованию:  нет

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

int main(void)
{
    /* Initialized data */
    static int in[9] = { 1,6,11,16,21,26,31,36,51 };

    /* Local variables */
    extern int ia83r_c(int *, int *, float *, float *, int *, int *,
                       float *, float *, int *, float *, int *);
    static int ierr, iter;
    static float a[8], h__;
    static int k, m, n, i;
    static float x[51], y[51];
    static int jn[9];
    static float ax[9], rab[18]  /* was [2][9] */;
    int i__1;
    float r__1;

#define rab_ref(a_1,a_2) rab[(a_2)*2 + a_1 - 3]

    m = 51;
    n = 7;
    h__ = 2.f / (m - 1);
    i__1 = m;
    for (k = 1; k <= i__1; ++k) {
        x[k - 1] = (k - 1) * h__ - 1.f;
        y[k - 1] = (r__1 = x[k - 1], abs(r__1));
/* L5: */
    }
    ia83r_c(&m, &n, x, y, in, &iter, a, ax, jn, rab, &ierr);

    printf("\n %5i \n",ierr);
    printf("\n %5i %5i %5i \n",n,m,iter);
    printf("\n %16.7e \n",rab_ref(1,1));
    for (i = 1; i <= 9; ++i) {
        printf("\n %16.7e \n",ax[i-1]);
    }
    for (i = 1; i <= 8; ++i) {
        printf("\n %16.7e \n",a[i-1]);
    }
    printf("\n jn = %3i %3i %3i %3i %3i %3i %3i %3i %3i \n",
               jn[0],jn[1],jn[2],jn[3],jn[4],jn[5],jn[6],jn[7],jn[8]);
    return 0;
} /* main */


Результаты:   
      
       ierr = 0;     iter = 6;
       rab(1,1) = 0.04587;
 
       a = (0.04587, 0., 2.8698, 0., - 4.18189, 0., 2.31212, 0.);
       ax = (- 1., -.880, -.56, -.2, 0., .2, .56, .88, 1.);

       jn = (1, 4, 12, 21, 26, 31, 40, 48, 51)