Текст подпрограммы и версий
ip03r_p.zip , ip03e_p.zip 		Тексты тестовых примеров
tip03r_p.zip , tip03e_p.zip

Подпрограмма:  IP03R (модуль IP03R_p)

Назначение

Полиномиальная интерполяция функции двух переменных, заданной на прямоугольной сетке, по модифицированной схеме Невилла.

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

Пусть заданы узлы неравномерной прямоугольной сетки: 

   X1A  =  { x1(1), x2(1), ..., xM(1) } ;    X2A  =  { x1(2), x2(2), ..., xN(2) }

 и значения  функции     f (x(1), x(2))

 в узлах этой сетки:     YA (I, J)  =  f (X1A(I), X2A(J)) ,

   I  =  1, 2, ..., M ;     J  =  1, 2, ..., N .

 Предполагается, что
    x1(1) < x2(1) < ... < xM(1)   и    x1(2) < x2(2) < ... < xN(2).

Подпрограмма IP03R вычисляет значение Y = f (X1, X2) в заданной точке плоскости с координатами (X1, X2) и оценку ошибки DY полученного значения. При этом по горизонтальному направлению используется интерполяционный полином степени MP, а по вертикальному направлению - интерполяционный полином степени NP. Вычисление соответствующих значений интерполяционных полиномов осуществляется по модифицированной схеме Невилла.

Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.

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

procedure IP03R(var X1A :Array of Real; var X2A :Array of Real;
                var YA :Array of Real; M :Integer; N :Integer;
                MP :Integer; NP :Integer; X1 :Real;
                X2 :Real; var Y :Real; var DY :Real;
                var YM :Array of Real; var YN :Array of Real;
                var C :Array of Real; var D :Array of Real;
                var IERR :Integer);


   Параметры
         X1A -
         X2A   		вещественные векторы длины M и N, содержащие узлы  { x1 (1), x2 (1), ..., xM (1) } и { x1 (2), x2 (2), ..., xN (2) } соответственно;
YA - вещественный двумерный массив размеров M на N, содержащий значения интерполируемой функции двух переменных в узлах заданной сетки;
M, N - длины векторов X1A и X2A соответственно, M ≥ MP + 1, N ≥ NP + 1 (тип: целый);
MP, NP - заданные степени интерполяционных полиномов, MP ≥ 1, NP ≥ 1 (тип: целый);
X1, X2 - координаты заданной точки на плоскости, в которой ищется значение интерполируемой функции (тип: вещественный);
Y - вещественная переменная, содержащая значение интерполируемой функции в заданной точке плоскости;
DY - вещественная переменная, содержащая оценку ошибки вычисленного значения интерполируемой функции;
YM, YN - вещественные векторы длины M и N соответственно, используемые в подпрограмме в качестве рабочих;
C, D - вещественные векторы длины  max (MP + 1, NP + 1), используемые в подпрограмме в качестве рабочих;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом
IERR=65 - когда по крайней мере два узла сетки по какому - либо направлению совпадают;
IERR=66 - когда число узлов сетки по какому - либо направлению меньше или равно степени соответствующего интерполяционного полинома. 
   Версии
IP03E - полиномиальная интерполяция функции двух переменных, заданной на прямоугольной неравномерной сетке, по модифицированной схеме Невилла в режиме расширенной (Extended) точности; при этом параметры X1A, X2A, YA, X1, X2, Y, DY, YM, YN, C, D должны иметь тип Extended. 
   Вызываемые подпрограммы
IP02R - вычисление значения интерполяционного полинома в заданной точке по модифицированной схеме Невилла. используется в подпрограмме IP03R;
IP02E - вычисление значения интерполяционного полинома в заданной точке по модифицированной схеме Невилла в режиме расширенной (Extended) точности; используется в подпрограмме IP03E. 
   Замечания по использованию
  Заданная точка на плоскости не обязательно должна лежать внутри заданной прямоугольной сетки. В этом случае подпрограмма IP03R и IP03E выполняют полиномиальную экстраполяцию двух переменных.

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

Unit TIP03R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IP03R_p;

function TIP03R: String;

implementation

function TIP03R: String;
var
M,N,I,J,MP,NP,IERR :Integer;
R,X1,X2,Y,DY :Real;
X1A :Array [0..14] of Real;
X2A :Array [0..9] of Real;
YA :Array [0..149] of Real;
YN :Array [0..9] of Real;
YM :Array [0..14] of Real;
C :Array [0..5] of Real;
D :Array [0..5] of Real;
label
_1,_2,_3;
begin
Result := '';  { результат функции }
M := 15;
N := 10;
R := 0.0;
for I:=1 to M do
 begin
  X1A[I-1] := R;
_1:
  R := R+0.1;
 end;
R := 0.0;
for I:=1 to N do
 begin
  X2A[I-1] := R;
_2:
  R := R+0.1;
 end;
for I:=1 to M do
 begin
  for J:=1 to N do
   begin
_3:
    YA[(I-1)+(J-1)*15] := Sin(X1A[I-1])*Sin(X2A[J-1]);
   end;
 end;
МР := 5;
NP := 4;
X1 := 0.55;
X2 := 0.65;
IP03R(X1A,X2A,YA,M,N,MP,NP,X1,X2,Y,DY,
     YM,YN,C,D,IERR);
Result := Result + Format(' %20.16f %20.16f ',[Y,DY]) + #$0D#$0A;
Result := Result + Format('%5d ',[IERR]) + #$0D#$0A;
UtRes('TIP03R',Result);  { вывод результатов в файл TIP03R.res }
exit;
end;

end.
     

 Результаты:
     
             Y = 0.316323 
          DY = 0.616967E-07 
       IERR = 0