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

Подпрограмма:  IIS3R (модуль IIS3R_p)

Назначение

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

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

Вычисляются коэффициенты кубического сплайна 

 S (x) = yi + c1 i ( x - xi ) + c2 i ( x - xi )2 + c3 i ( x - xi )3 ,
 x ∈ [xi, xi + 1]  ,    i = 1, 2, ..., N - 1  ,

соответствующего сетке  {x i} и интерполирующего заданные значения табличной функции  yi = y (x i) ,  i = 1, 2, ..., N в узлах сетки.

Коэффициенты  c1 i , c2 i , c3 i ,  i = 1, 2, ..., N - 1 определяются значениями  S''i = S'' (x i) ,  i = 1, 2, ..., N, для нахождения которых решается система линейных алгебраических уравнений, полученная из условия непрерывности функции  S' (x) в узлах сетки и из краевых условий вида: 

          2 * S''1 + b1 * S''2 = b2
          b3 * S''N - 1 + 2 * S''N = b4  ,

где  b1 , b2 , b3 , b4 - заданные параметры. Если значения производных  y''1 = y'' (x1) и  y''N = y'' (xN) интерполируемой функции известны, то надо положить: 

          b1 = 0 ,  b2 = 2 * y''1
          b3 = 0 ,  b4 = 2 * y''N

Если известны значения  y'1 = y' (x1) и  y'N = y' (xN), то надо положить: 

 b1 = 1 ,   b2 = (6/(x2 - x1)) * ( (y2 - y1)/(x2 - x1) - y'1 )  ,
 b3 = 1 ,   b4 = (6/(xN - xN - 1)) * ( y'N - (yN - yN - 1)/(xN - xN - 1) )

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

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

procedure IIS3R(var X :Array of Real; var Y :Array of Real;
                NX :Integer; var BPAR :Array of Real;
                var C :Array of Real; var IERR :Integer);


   Параметры
X - вещественный вектоp длины  NX, содержащий заданные значения узлов сетки и упорядоченный так, что  X (1) < X (2) < ... < X (NX);
Y - вещественный вектоp длины  NX, содержащий заданные в узлах сетки значения интерполируемой функции;
NX - заданное число узлов сетки,  NX ≥ 2 (тип: целый);
BPAR - вещественный вектоp длины 4, содержащий параметры краевых условий,  BPAR (I) = bI,  I = 1, 2, 3, 4;
C - вещественный двумерный массив размера  3 * (NX - 1), содержащий вычисленные коэффициенты интерполяционного кубического сплайна;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
IERR=65 - когда заданное число узлов сетки меньше 2;
IERR=66 - когда элементы вектоpа  X не упорядочены по возрастанию. 
   Версии:  нет 
   Вызываемые подпрограммы
UTI I10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы I IS3R. 
   Замечания по использованию:  нет

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

Unit tiis3r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IIS3R_p;

function tiis3r: String;

implementation

function tiis3r: String;
var
NX,I,J,IERR :Integer;
Y :Array [0..3] of Real;
BPAR :Array [0..3] of Real;
C :Array [0..8] of Real;
const
X :Array [0..3] of Real = ( 0.0,0.33,1.0,2.0 );
label
_5,_10;
begin
Result := '';  { результат функции }
NX := 4;
for I:=1 to NX do
 begin
  Y[I-1] := ((X[I-1]-2.0)*X[I-1]+3.0)*X[I-1]+4.0;
_5:
 end;
BPAR[0] := 1.0;
BPAR[1] := 6.0/(X[1]-X[0])*((Y[1]-Y[0])/(X[1]-X[0])-3.0);
BPAR[2] := 1.0;
BPAR[3] := 6.0/(X[3]-X[2])*(7.0-(Y[3]-Y[2])/(X[3]-X[2]));
IIS3R(X,Y,NX,BPAR,C,IERR);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
for I:=1 to 3 do
 begin
_10:
  Result := Result + #$0D#$0A;
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[C[(I-1)+(J-1)*3]]) + #$0D#$0A;
   end;
  Result := Result + #$0D#$0A;
 end;
UtRes('tiis3r',Result);  { вывод результатов в файл tiis3r.res }
exit;
end;

end.


Результаты:

      IERR = 0
 
                 |  3.0000  -2.0000  1.0000  |
      C  =   |  2.0067  -1.0100  1.0000  |
                |  2.0000   1.0000  1.0000  |