Текст подпрограммы и версий
is01r_p.zip , is01e_p.zip 		Тексты тестовых примеров
tis01r_p.zip , tis01e_p.zip

Подпрограмма:  IS01R (модуль IS01R_p)

Назначение

Построение одномерного сглаживающего кубического сплайна.

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

Пусть на сетке x1 < x2 < ... < xN заданы измеренные значения f1, f2, ..., fN некоторой гладкой функции. Требуется найти дважды непрерывно дифференцируемую функцию  s (x), которая минимизирует функционал 

        xN
        ∫   ( g''(x) )2 dx
       xi
  среди всех функций g (x) таких, что
         N 
        ∑    [ ( g(xi) - fi ) / ( dfi ) ]2  ≤  SM ,
        i=1

где  dfi > 0 и SM ≥ 0 - заданные числа.

Числа  dfi характеризуют точность задания дискретных данных  fi, а SM - требуемую меpу их сглаживания.

Решением задачи является сглаживающий кубический сплайн, который на каждом интервале [xi, xi + 1] представляется в виде 

   s(x)  =  c3 i h3 + c2 i h2 + c1 i h + yi  ,    h = x - xi ,

где  c3 i , c2 i , c1 i ,  i = 1, 2, ..., N - 1 суть коэффициенты, а  yi - значения сглаживающего сплайна в узлах  xi ,  i = 1, 2, ..., N.

Reinsch C.H. Smoothing by Spline Functions. Numerishe Mathematik, 1967, v.10, 177 - 183.

Морозов B.A. O задаче дифференцирования и некоторых алгоритмах приближения экспериментальной информации, в сб. "Вычислительные методы и программирование", 1970, вып.ХIV, Изд - во МГУ, 46 - 62.

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

procedure IS01R(var X :Array of Real; var F :Array of Real;
                var DF :Array of Real; NX :Integer; SM :Real;
                var Y :Array of Real; var C :Array of Real;
                var RAB :Array of Real; var IERR :Integer);


   Параметры
X - вещественный массив длины NX, в котоpом задаются значения узлов сетки такие, что X (I) < X (I + 1),  I = 1, ..., NX - 1;
F - вещественный массив длины NX, в котоpом задаются значения сглаживаемой дискретной функции;
DF - вещественный массив длины NX, в котоpом задаются значения весов DF (I) > 0 ,  I = 1, ..., NX;
NX - заданное число узлов, NX ≥ 2 (тип: целый);
SM - неотрицательное число, задающее меpу сглаживания (тип: вещественный);
Y - вещественный массив длины NX значений сглаживающего сплайна в узлах сетки;
C - двумерный вещественный массив размера 3* (NX - 1) значений коэффициентов сглаживающего сплайна;
RAB - двумерный рабочий массив размера 7* (NX + 2) (тип: вещественный);
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы:
IERR=65 - если NX < 2;
IERR=66 - если массив X не упорядочен строго по возрастанию. 
   Версии
IS01E - построение одномерного сглаживающего кубического сплайна с расширенной (Extended) точностью. 
   Вызываемые подпрограммы
UTIS10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы IS01R.
UTIS11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы IS01E. 
   Замечания по использованию:  нет

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

Unit TIS01R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IS01R_p;

function TIS01R: String;

implementation

function TIS01R: String;
var
NX,_i,K,I,IERR :Integer;
SM :Real;
Y :Array [0..4] of Real;
C :Array [0..11] of Real;
RАВ :Array [0..48] of Real;
const
X :Array [0..4] of Real = ( 1.0,2.0,3.0,4.0,5.0 );
F :Array [0..4] of Real = ( 1.0,4.0,20.0,16.0,25.0 );
DF :Array [0..4] of Real = ( 0.2,0.2,10.0,0.2,0.2 );
begin
Result := '';  { результат функции }
NX := 5;
SM := NX-1.0-Sqrt((NX+NX));
IS01R(X,F,DF,NX,SM,Y,C,RAB,IERR);
Result := Result + Format('%s',[' IERR']);
Result := Result + Format('%2d ',[IERR]);
Result := Result + Format('%s',[#$0D#$0A + ' Y']);
Result := Result + #$0D#$0A;
for _i:=0 to 4 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to 3 do
 begin
  for K:=1 to 4 do
   begin
    Result := Result + Format(' C%20.16f ',
 [C[(I-1)+(K-1)*3]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
UtRes('TIS01R',Result);  { вывод результатов в файл TIS01R.res }
exit;
end;

end.

Результаты:

       IERR  =  0

       Y(1) = .999 ; Y(2) = 4.003 ; Y(3) = 10.85 ; Y(4) = 16.003 ; Y(5) = 24.999

       C(1, 1) =  1.785   C(1, 2) =  5.442   C(1, 3) =  6.000   C(1, 4) =  6.558
       C(2, 1) =  0.000   C(2, 2) =  3.657   C(2, 3) = -3.099   C(2, 4) =  3.657
       C(3, 1) =  1.219   C(3, 2) = -2.252   C(3, 3) =  2.252   C(3, 4) = -1.219