БЧА НИВЦ МГУ. IIB9R_P. Сплайн-интерполяция

Текст подпрограммы и версий
iib9r_p.zip , iib9e_p.zip

Тексты тестовых примеров
tiib9r_p.zip , tiib9e_p.zip

Подпрограмма: IIB9R (модуль IIB9R_p)

Назначение

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

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

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

   X1A  =  { x₁⁽¹⁾, x₂⁽¹⁾, ..., x_M⁽¹⁾ } ;   X2A  =  { x₁⁽²⁾, x₂⁽²⁾, ..., x_N⁽²⁾ } ,

 значения интерполируемой функции    f (x⁽¹⁾, x⁽²⁾)

 и ее второй частной производной   ∂⁽²⁾ f (x⁽¹⁾, x⁽²⁾) / ∂x^{(2) 2}
 в узлах этой сетки:
 YA(I, J)  =  f (X1A(I), X2A(J)) ;   Y2A(I, J)  =  ∂² f (X1A(I), X2A(J)) / ∂x^{(2) 2} ,

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

а также следующие значения первой частной производной функции f по первой переменной:

   YP1  =  ∂ f (X1A(1), X2)) / ∂x⁽¹⁾ ,   YPM  =  ∂ f (X1A(M), X2) / ∂x⁽¹⁾ .

Здесь X2 - вторая координата заданной точки (X1, X2) внутри сеточного прямоугольника, в которой ищется значение интерполируемой функции: Y = f (X1, X2). Используется метод бикубических сплайнов.

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

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

procedure IIB9R(var X1A :Array of Real; var X2A :Array of Real;
                var YA :Array of Real; var Y2A :Array of Real;
                M :Integer; N :Integer; X1 :Real;
                X2 :Real; YP1 :Real; YPM :Real;
                var Y :Real; var YT :Array of Real;
                var Y2T :Array of Real; var YYT :Array of Real;
                var U :Array of Real);

   Параметры

X1A - X2A	вещественные векторы длины M и N, содержащие узлы {x₁ ⁽¹⁾, x₂ ⁽¹⁾, ...,x_M ⁽¹⁾} и {x₁ ⁽²⁾, x₂ ⁽²⁾, ...,x_N ⁽²⁾} соответственно;
YA -	вещественный двумерный массив размеров M на N, содержащий значения интерполируемой функции f (x⁽¹⁾, x⁽²⁾) двух переменных в узлах заданной сетки;
Y2A -	вещественный двумерный массив размеров M и N, содержащий значения второй частной производной функции f по второй переменной в узлах заданной сетки;
M, N -	длины векторов X1A и X2A соответственно (тип: целый);
X1, X2 -	координаты заданной точки внутри сеточного прямоугольника, в которой ищется значение интерполируемой функции f (тип: вещественный);
YP1 - YPM	заданные значения первой частной производной функции f по первой переменной в точках (X1A (1), X2) и (X1A (M), X2) соответственно (тип: вещественный);
Y -	вещественная переменная, содержащая вычисленное значение интерполируемой функции f в точке (X1, X2);
YT -	вещественный вектор длины N, используемый в подпрограмме в качестве рабочего;
Y2T -	вещественный вектор длины max (M, N), используемый в подпрограмме в качестве рабочего;
YYT, U -	вещественные векторы длины M, используемые в подпрограмме в качестве рабочих.

   Версии

I IB9E -

выполняет те же вычисления, что и подпрограмма I IB9R, однако в режиме расширенной (Extended) точности. При этом все формальные параметры, кроме M и N, должны иметь тип Extended.

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

IDS8R -
IDS8E

вычисление вторых производных таблично - заданной функции в узлах неравномерной сетки методом кубических или натуральных кубических сплайнов при заданных первых производных функции в концевых узлах сетки в режиме одинарной и расширенной (Extended) точности соответственно; используются в подпрограммах I IB9R и I IB9E соответственно.

I IS8R -
I IS8E

вычисление значения в заданной точке интерполируемой табличной функции, определенной в узлах неравномерной сетки, упорядоченной по возрастанию, по известным значениям ее второй производной в узлах этой сетки методом кубических сплайнов в режиме одинарной и расширенной (Extended) точности соответственно; используются в подпрограммах I IB9R и I IB9E соответственно.

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

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

Unit TIIB9R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IIB9R_p;

function TIIB9R: String;

implementation

function TIIB9R: String;
var
M,N,I,J :Integer;
R,X1,X2,YP1,YPM,Y1T,Y :Real;
X1A :Array [0..4] of Real;
X2A :Array [0..3] of Real;
YA :Array [0..19] of Real;
Y2A :Array [0..19] of Real;
YT :Array [0..3] of Real;
Y2T :Array [0..4] of Real;
YYT :Array [0..4] of Real;
U :Array [0..4] of Real;
label
_1,_2,_4;
begin
Result := '';  { результат функции }
M := 5;
N := 4;
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
    YA[(I-1)+(J-1)*5] := Sin(X1A[I-1])*Sin(X2A[J-1]);
_4:
    Y2A[(I-1)+(J-1)*5] := -YA[(I-1)+(J-1)*5];
   end;
 end;
X1 := 0.55;
X2 := 0.65;
YP1 := Cos(X1A[0])*Sin(X2);
YРМ := Cos(X1A[M-1])*Sin(X2);
IIB9R(X1A,X2A,YA,Y2A,M,N,X1,X2,YP1,YPM,Y,
     YT,Y2T,YYT,U);
Y1T := Sin(X1)*Sin(X2);
Result := Result + Format(' %20.16f %20.16f ',[Y,Y1T]) + #$0D#$0A;
UtRes('TIIB9R',Result);  { вывод результатов в файл TIIB9R.res }
exit;
end;

end.
      

Результат:   Y  =  0.316370