Текст подпрограммы и версий ( Фортран )
ip03r.zip , ip03d.zip
Тексты тестовых примеров ( Фортран )
tip03r.zip , tip03d.zip
Текст подпрограммы и версий ( Си )
ip03r_c.zip , ip03d_c.zip
Тексты тестовых примеров ( Си )
tip03r_c.zip , tip03d_c.zip
Текст подпрограммы и версий ( Паскаль )
ip03r_p.zip , ip03e_p.zip
Тексты тестовых примеров ( Паскаль )
tip03r_p.zip , tip03e_p.zip

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

Назначение

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

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

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

   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.

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

    SUBROUTINE  IP03R (X1A, X2A, YA, M, N, MP, NP, X1, X2, Y, DY,
                                          YM, YN, C, D, IERR) 

Параметры

         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 - когда число узлов сетки по какому - либо направлению меньше или равно степени соответствующего интерполяционного полинома.

Версии

IP03D - полиномиальная интерполяция функции двух переменных, заданной на прямоугольной неравномерной сетке, по модифицированной схеме Невилла в режиме удвоенной точности; при этом параметры X1A, X2A, YA, X1, X2, Y, DY, YM, YN, C, D должны иметь тип DOUBLE PRECISION.

Вызываемые подпрограммы

IP02R - вычисление значения интерполяционного полинома в заданной точке по модифицированной схеме Невилла; используется в подпрограмме IP03R;
IP02D - вычисление значения интерполяционного полинома в заданной точке по модифицированной схеме Невилла в режиме удвоенной точности; используется в подпрограмме IP03D.

Замечания по использованию

  Заданная точка на плоскости не обязательно должна лежать внутри заданной прямоугольной сетки. В этом случае подпрограмма IP03R и IP03D выполняют полиномиальную экстраполяцию двух переменных.

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

       DIMENSION  X1A(15), X2A(10), YA(15, 10), YM(15),
      *                     YN(10), C(6), D(6)
       M = 15 
       N = 10
       R = 0.0 
       DO 1  I = 1, M 
       X1A(I) = R 
    1 R = R + 0.1 
       R = 0.0 
       DO 2  I = 1, N 
       X2A(I) = R 
    2 R = R + 0.1 
       DO 3  I = 1, M 
       DO 3  J = 1, N 
    3 YA(I, J) = SIN(X1A(I))*SIN(X2A(J)) 
       MP = 5 
       NP = 4 
       X1 = 0.55 
       X2 = 0.65 
       CALL  IP03R (X1A, X2A, YA, M, N, MP, NP, X1, X2, Y, DY, YM,
      *                       YN, C, D, IERR)
     
 Результаты:
     
             Y = 0.316323 
          DY = 0.616967E-07 
       IERR = 0