Текст подпрограммы и версий
ia83r_p.zip , ia83e_p.zip 		Тексты тестовых примеров
tia83r_p.zip , tia83e_p.zip

Подпрограмма:  IA83R (модуль IA83R_p)

Назначение

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

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

Пусть в узлах сетки  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.

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

procedure IA83R(M :Integer; N :Integer; var X :Array of Real;
                var Y :Array of Real; var IN_ :Array of Integer;
                var ITER :Integer; var A :Array of Real;
                var AX :Array of Real; var JN :Array of Integer;
                var RAB :Array of Real; var IERR :Integer);


   Параметры
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 итераций. 
   Версии
IA83E - построение алгебраического полинома наилучшего равномерного приближения для таблично заданной функции в режиме вычислений с расширенной (Extended) точностью. 
   Вызываемые подпрограммы
UTIA10 - подпрограмма выдачи диагностических сообщений при работе IA83R.
UTIA11 - подпрограмма выдачи диагностических сообщений при работе IA83E. 
   Замечания по использованию:  нет

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

Unit TIA83R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IA83R_p;

function TIA83R: String;

implementation

function TIA83R: String;
var
M,N,K,I,ITER,IERR :Integer;
H :Real;
A :Array [0..7] of Real;
АХ :Array [0..8] of Real;
RАВ :Array [0..17] of Real;
JN :Array [0..8] of Integer;
X :Array [0..50] of Real;
Y :Array [0..50] of Real;
const
IN_ :Array [0..8] of Integer = ( 1,6,11,16,21,26,31,36,51 );
label
_5;
begin
Result := '';  { результат функции }
M := 51;
N := 7;
H := 2.0/(M-1);
for K:=1 to M do
 begin
  X[K-1] := -1.0+(K-1)*H;
  Y[K-1] := Abs(X[K-1]);
_5:
 end;
IA83R(M,N,X,Y,IN_,ITER,A,AX,JN,RAB,IERR);
Result := Result + Format('%s',
 [' ТЕСТОВЫЙ ПРИМЕР ДЛЯ IA83R' + #$0D#$0A]) + #$0D#$0A; 
Result := Result + Format('%s',[' PEЗYЛЬTATЫ:' + #$0D#$0A + ' IERR=']);
Result := Result + Format('%2d ',[IERR]);
Result := Result + Format('%s',['   N=']);
Result := Result + Format('%2d ',[N]);
Result := Result + Format('%s',['   M=']);
Result := Result + Format('%3d ',[M]);
Result := Result + Format('%s',['   ITER=']);
Result := Result + Format('%3d ',[ITER]);
Result := Result + Format('%s',['   MINMAX=']);
Result := Result + Format('%20.16f ',[RAB[0]]) + #$0D#$0A;
for I:=1 to 8 do
 begin
  Result := Result + Format('       A(%1d )=%20.16f ',[I,A[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
  for I:=1 to 9 do
   begin
    Result := Result + Format(' АХ(%1d )=%20.16f ',[I,AX[I-1]]) + #$0D#$0A;
   end;
Result := Result + #$0D#$0A;
    for I:=1 to 9 do
     begin
      Result := Result + Format(' JN(%1d )=%5d ',[I,JN[I-1]]) + #$0D#$0A;
     end;
Result := Result + #$0D#$0A;
UtRes('TIA83R',Result);  { вывод результатов в файл TIA83R.res }
exit;
end;

end.

Результаты:   
      
       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)