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

Подпрограмма:  IS02R (модуль IS02R_p)

Назначение

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

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

Пусть в узлах неравномерной сетки X1 < X2 < ... < XN приближенно заданы дискретные значения Uk,  k = 1, 2, ..., N. Обозначим через D множество дважды непрерывно дифференцируемых функций  g (x), удовлетворяющих условиям периодичности: 

   g(u) (X1)  =  g(u) (XN)  ,   u = 0, 1, 2 .

Требуется найти сглаженные значения y (Xk),  y" (Xk) периодического кубического сплайна  y (X), минимизирующего функционал 

      N
     ∑   pk( g(xk) - uk )2  +
    k=1
                                     XN
                                +   ∫   ( g"(x) )2 dx ,   g ∈ D
                                   X1

Здесь  pk > 0,  k = 1, 2, ..., N, суть заданные веса, определяемые точностью задания  uk и требуемой мерой их сглаживания. В частности, при 

 pk = σk2 / α   ,    σk2 = M ( uk - M uk )2 ,

где M - символ математического ожидания, получаем регуляризованный периодический кубический сплайн  y (x), соответствующий параметру регуляризации α > 0.

Для контроля за степенью сглаживания вычисляются такие характеристики: 

             XN
   E1  =   ∫   [ y"(x) ]2 dx  ,
             X1
                  N
   E2  =   (  ∑  pk ( yk - uk )2 ) / ( p1 + ... + pN )
                 k=1

Сплайн  y (x) на отрезке [Xk, Xk + 1] может быть восстановлен по формуле 

                                (xk+1-x)3                              (x-xk)3
   y(x)  =  y"(xk) * -------------   +   y"(xk+1)  *  -------------  +
                                 6*δxk                                  6*δxk

                              y"(xk)*(δxk)2                   xk+1-x
    +  ( y(xk)  -   -----------------------  )  *  -----------------  +
                                       6                               δxk

                               y"(xk+1)*(δxk)2               x-xk
    +  ( y(xk+1)  -  ----------------------  )  * ---------------   ,
                                          6                             xk

 где     Δxk = xk + 1 - xk .

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

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

procedure IS02R(N :Integer; var X :Array of Real;
                var U :Array of Real; var P :Array of Real;
                var Y :Array of Real; var Y2 :Array of Real;
                var DX :Array of Real; var E :Array of Real;
                var R :Array of Real; var IERR :Integer);


   Параметры
N - заданное число узлов сетки, N ≥ 3 (тип: целый);
X - вещественный вектоp длины N, содержащий заданные значения узлов сетки и упорядоченный так, что X (1) < X (2) < ...< X (N);
U - вещественный вектоp длины N, содержащий заданные значения сглаживаемой функции;
P - вещественный вектоp длины N, содержащий заданные веса, P (I) > 0, I = 1, 2, ..., N;
Y - вещественный вектоp длины N, содержащий вычисленные в узлах сетки значения сглаживающего периодического кубического сплайна;
Y2 - вещественный вектоp длины N, содержащий вычисленные в узлах сетки значения второй производной сглаживающего периодического кубического сплайна;
DX - вещественный вектоp длины N - 1, содержащий вычисленные значения шагов сетки  DX (K) = X (K + 1) - X (K),  K = 1, 2, ..., N - 1;
E - вещественный вектоp длины 2,  E (1) = E1,  E (2) = E2;
R - вещественный двумерный рабочий массив размера 6 на N;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
IERR=65 - если N < 3;
IERR=66 - если элементы вектоpа X не упорядочены строго по возрастанию;
IERR=67 - если хотя бы один элемент вектоpа P меньше или pавен нулю. 
   Версии:  нет 
   Вызываемые подпрограммы
UTI I10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы IS02R. 
   Замечания по использованию
  Подпрограмма IS02R может быть использована для сглаживания приближенно заданных на неравномерной сетке t1 < t2 < ...< tN значений  uk = u (tk),  vk = v (tk),  k = 1, 2, ..., N параметрической циклической функции  u = u (t),  v = v (t). Для этого достаточно дважды обратиться к подпрограмме IS02R.

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

Unit tis02r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IS02R_p;

function tis02r: String;

implementation

function tis02r: String;
var
N,I,_i,IERR :Integer;
PI,PI2,H :Real;
X :Array [0..8] of Real;
U :Array [0..8] of Real;
P :Array [0..8] of Real;
Y :Array [0..8] of Real;
Y2 :Array [0..8] of Real;
DX :Array [0..7] of Real;
E :Array [0..1] of Real;
R :Array [0..53] of Real;
label
_1;
begin
Result := '';  { результат функции }
PI := 3.141592;
PI2 := 2*PI;
N := 9;
H := PI2/(N-1);
X[0] := 0.0;
U[0] := 0.0;
P[0] := 1.0;
for I:=2 to N do
 begin
  X[I-1] := X[I-2]+H;
  U[I-1] := Sin(X[I-1]);
  P[I-1] := 1.0;
_1:
 end;
IS02R(N,X,U,P,Y,Y2,DX,E,R,IERR);
Result := Result + Format('%5d ',[IERR]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 8 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:=0 to 8 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;
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[DX[_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 1 do
 begin
  Result := Result + Format('%20.16f ',[E[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('tis02r',Result);  { вывод результатов в файл tis02r.res }
exit;
end;

end.


Результаты:

       IERR  =  0

       Y(1)  =   0.002020      Y2(1)  =   0.14309
       Y(2)  =   0.392894      Y2(2)  =  -0.410578
       Y(3)  =   0.555974      Y2(3)  =  -0.588419
       Y(4)  =   0.391777      Y2(4)  =  -0.417586
       Y(5)  =  -0.004407      Y2(5)  =   0.001235
       Y(6)  =  -0.399604      Y2(6)  =   0.423318
       Y(7)  =  -0.558361      Y2(7)  =   0.603963
       Y(8)  =  -0.380293      Y2(8)  =   0.437473
       Y(9)  =   0.002020      Y2(9)  =   0.014309

       DX(I)  =  0.785398 ,      I  =  1, 2, ..., 8
       E(1)  =  1.009781
       E(2)  =  0.087973