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

Подпрограмма:  IDS4R (модуль IDS4R_p)

Назначение

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

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

Пусть заданы пары значений (xi, yi). Считается, что  xi = x (ti),  yi = y (ti) ,  i = 1, 2, ... N, а величины ti = ti - 1 + di,  i =  2, 3, ... N,  t1 = 0. Шаги по параметру  di могут определятся из следующих формул: 

     а)  di  =  [ (Δxi)2  +  (Δyi)2 ]1/2

     б)  di  =  ( Δxi )2  +  ( Δyi )2

     в)  di  =  | Δxi |  +  | Δyi | ,   i = 2, ..., N-1 ,

где  Δ xi = xi + 1 - xi ,  Δ yi = yi + 1 - yi .

Далее по парам значений (ti, xi),  (ti, yi) ,  i = 1, ..., N и известным значениям вторых производных в концевых точках  x''1,  x''N и  y''1,  y''N методом кубических сплайнов [1] вычисляются значения вторых производных  x''i и  y''i,   i = 2, ...,N - 1.

Интерполирующие функции  x = x (t) и  y = y (t) на каждом из отрезков [ti, ti + 1] представляют из себя кубические полиномы, построенные по значениям (ti, xi),  (ti, x''i),  (ti + 1,xi + 1),  (ti + 1,x''i + 1) и  (ti, yi),  (ti, y''i),  (ti + 1, yi + 1),  (ti + 1, y''i + 1) как интерполяционные полиномы Эрмита [2]. Tак построенные функции  x (t) и  y (t) являются дважды непрерывно дифференцируемыми на всем отрезке [t1, tN] по параметру  t.

1. 

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

2.  И.С.Березин, Н.П.Жидков, Методы вычислений, т.1, Физматгиз, 1962.

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

procedure IDS4R(N :Integer; NPAR :Integer; var T :Array of Real;
                var X :Array of Real; var X2 :Array of Real;
                var Y :Array of Real; var Y2 :Array of Real;
                var F :Array of Real; var G :Array of Real);


   Параметры
N - число узлов сетки, N ≥ 3 (тип: целый);
NPAR - одно из целых чисел 1, 2, 3, характеризующих способ параметризации; при NPAR = 1 величины  di вычисляются по формуле а), при NPAR = 2 - по формуле б), при NPAR = 3 - по формуле в);
T - одномерный массив значений  ti размерности N (тип: вещественный);
X - одномерный массив значений  xi размерности N (тип: вещественный);
X2 - одномерный вещественный массив вычисляемых значений вторых производных функции  x (t) в узлах сетки  ti; на входе  X2 (1) = x''1,  X2 (N) = x''N; длина масива X2 pавна N;
Y - одномерный массив значений  yi размерности N (тип: вещественный);
Y2 - одномерный вещественный массив вычисляемых значений вторых производных функции  y (t) в узлах сетки  ti; на входе Y2 (1) = y''1,  Y2 (N) = y''N; длина массива Y2 pавна N;
F, G - одномерные рабочие массивы размерности N (тип: вещественный). 
   Версии:  нет 
   Вызываемые подпрограммы
IDS0R - вычисление значений вторых производных вещественной табличной функции от одного переменного в узлах неравномерной сетки методом интерполирующих кубических сплайнов при известных значениях вторых производных в концевых точках. 
   Замечания по использованию:  нет

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

На верхней части единичной полуокружности были взяты десять равноотстоящих точек. Вычислялись значения вторых производных  x''i и  y''i ,  i = 2, ..., 9 . 

Unit tids4r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IDS4R_p;

function tids4r: String;

implementation

function tids4r: String;
var
I,_i :Integer;
T :Array [0..9] of Real;
X :Array [0..9] of Real;
X2 :Array [0..9] of Real;
Y :Array [0..9] of Real;
Y2 :Array [0..9] of Real;
F :Array [0..9] of Real;
G :Array [0..9] of Real;
label
_5;
begin
Result := '';  { результат функции }
for I:=1 to 10 do
 begin
  X[I-1] := -1+(I-1)*(2.0/9.0);
_5:
  Y[I-1] := Sqrt(1-X[I-1]*X[I-1]);
 end;
X2[0] := 1.0;
X2[9] := -1;
Y2[0] := 0.0;
Y2[9] := 0.0;
IDS4R(10,1,T,X,X2,Y,Y2,F,G);
Result := Result + #$0D#$0A;
for _i:=0 to 9 do
 begin
  Result := Result + Format('%20.16f ',[X2[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 9 do
 begin
  Result := Result + Format('%20.16f ',[Y2[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('tids4r',Result);  { вывод результатов в файл tids4r.res }
exit;
end;

end.


Результат:

       X2  =  ( 1.00,   0.83,   0.55,   0.34,   0.11,  -0.11,  -0.34, 
                    -0.55,  -0.83,  -1.00 )
       Y2  =  ( 0.00,  -0.70,  -0.83,  -0.96,  -1.00,  -1.00,  -0.96, 
                    -0.83,  -0.70,   0.00 )