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

Подпрограмма:  IIS6R (модуль IIS6R_p)

Назначение

Интерполяция на заданном элементарном прямоугольнике натуральным бикубическим сплайном табличной функции от двух переменных.

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

Hа заданном элементарном прямоугольнике  R = [xlx, xlx + 1] * [yly, yly + 1] вычисляются коэффициенты  CRm n натурального бикубического сплайна 

                          4       4            
        S (x, y)  =  ∑     ∑   CRm n (x - xlx)m -1 (y - yly)n - 1 ,
                       m = 1  n = 1
               (x, y) ∈ R  ,   1 ≤ lx ≤ NX  ,   1 ≤ ly ≤ NY  ,

соответствующего двумерной сетке {x i, yj} и интерполирующего заданные значения  fi j = f (x i, yj),  i = 1, 2, ..., NX,   j = 1, 2, ..., NY, в узлах сетки. Для нахождения 16 - ти коэффициентов  CRm n,  m, n = 1, 2, 3, 4, решается система линейных алгебраических уравнений, полученная по значениям бикубического сплайна в 16 - ти точках заданного прямоугольника.

Эти значения вычисляет библиотечная подпрограмма I IS5R.

Дж. Алберг, Э.Нильсон, Дж.Уолш. Теория сплайнов и ее приложения, M., Мир, 1972.

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

procedure IIS6R(var F :Array of Real; var X :Array of Real;
                NX :Integer; var Y :Array of Real;
                NY :Integer; LX :Integer; LY :Integer;
                var C :Array of Real; var WK :Array of Real;
                var IERR :Integer);


   Параметры
F - вещественный двумерный массив размера  NX * NY, содержащий заданные в узлах двумерной сетки значения интерполируемой функции;
X - вещественный вектоp длины  NX, содержащий заданные значения узлов одномерной сетки по оси  x и упорядоченный так, что  X (1) < X (2) < ... < X (NX);
NX - заданное число узлов одномерной сетки по оси  x;   NX ≥ 2 (тип: целый);
Y - вещественный вектоp длины  NY, содержащий заданные значения узлов одномерной сетки по оси  y и упорядоченный так, что  Y (1) < Y (2) < ... < Y (NY);
NY - заданное число узлов одномерной сетки по оси  y;   NY ≥ 2 (тип: целый);
LX - заданное число, определяющее элементарный прямоугольник  R по оси  x,  1 ≤ LX < NX (тип: целый);
LY - заданное число, определяющее элементарный прямоугольник  R по оси  y,  1 ≤ LY < NY (тип: целый);
C - вещественный двумерный массив размера 4 * 4, содержащий вычисленные на заданном элементарном прямоугольнике  R = [x lx, x lx + 1] * [yly, yly + 1] коэффициенты интерполяционного бикубического сплайна;
WK - вещественный вектоp длины  4 * MAX (4, NY) + MAX ((NX - 1) * 3, (NY - 1) * 3 + 4), используемый подпрограммой как рабочий;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
IERR=65 - когда заданное число узлов сетки по  x меньше 2;
IERR=66 - когда заданное число узлов сетки по  y меньше 2;
IERR=67 - когда элементы вектоpа  X не упорядочены по возрастанию;
IERR=68 - когда элементы вектоpа  Y не упорядочены по возрастанию;
IERR=69 - когда  LX ≥ NX или  LX < 1;
IERR=70 - когда  LY ≥ NY или  LY < 1. 
   Версии:  нет 
   Вызываемые подпрограммы
I IS5R - вычисление значений интерполяционного натурального бикубического сплайна на заданном множестве точек.
I IS4R - вычисление значений интерполяционного кубического сплайна в заданных точках.
I IS3R - интерполяция кубическими сплайнами табличной функции от одной переменной на неравномерной сетке при известных краевых условиях на вторые производные.
UTI I10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы I IS6R. 
   Замечания по использованию:  нет

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

Unit tiis6r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IIS6R_p;

function tiis6r: String;

implementation

function tiis6r: String;
var
NX,NY,LX,LY,I,J,IERR :Integer;
C :Array [0..15] of Real;
WK :Array [0..25] of Real;
const
X :Array [0..2] of Real = ( 1.0,3.0,5.0 );
Y :Array [0..2] of Real = ( 1.0,3.0,5.0 );
F :Array [0..8] of Real = ( 2.0,4.0,3.0,0.5,1.0,0.75,1.0,2.0,1.5 );
label
_5;
begin
Result := '';  { результат функции }
NX := 3;
NY := 3;
LX := 1;
LY := 1;
IIS6R(F,X,NX,Y,NY,LX,LY,C,WK,IERR);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
for I:=1 to 4 do
 begin
_5:
  Result := Result + #$0D#$0A;
  for J:=1 to 4 do
   begin
    Result := Result + Format('%20.16f ',[C[(I-1)+(J-1)*4]]) + #$0D#$0A;
   end;
  Result := Result + #$0D#$0A;
 end;
UtRes('tiis6r',Result);  { вывод результатов в файл tiis6r.res }
exit;
end;

end.


Результаты:

       IERR  =  0

                |  2.000000  -1.000000   0.000000    0.062500  |
       C  =  |  1.375000  -0.687500   0.000000    0.042969  |
                |  0.000000   0.000000   0.000000    0.000000  |
                | -0.093750   0.046875   0.000000  -0.002930  |

   т.е. для  (x, y) ∈ [1., 3.] * [1., 3.]

                          4       4            
        S (x, y)  =  ∈     ∑   Cm n (x - 1.)m -1 (y - 1.)n - 1 
                       m = 1  n = 1