Текст подпрограммы и версий ( Фортран )
is02r.zip 		Тексты тестовых примеров ( Фортран )
tis02r.zip
Текст подпрограммы и версий ( Си )
is02r_c.zip 		Тексты тестовых примеров ( Си )
tis02r_c.zip
Текст подпрограммы и версий ( Паскаль )
is02r_p.zip 		Тексты тестовых примеров ( Паскаль )
tis02r_p.zip

Подпрограмма:  IS02R

Назначение

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

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

Пусть в узлах неравномерной сетки 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.

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

    SUBROUTINE  IS02R (N, X, U, P, Y, Y2, DX, E, R, IERR) 
   Параметры
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.

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

       DIMENSION  X(9), U(9), P(9), Y(9), Y2(9), DX(8), E(2), R(6, 9)
       PI = 3.141592
       PI2 = 2*PI
       N = 9
       H = PI2/(N - 1)
       X(1) = 0.
       U(1) = 0.
       P(1) = 1.
       DO 1  I = 2, N
       X(I) = X(I-1) + H
       U(I) = SIN( X(I) )
       P(I) = 1.
    1 CONTINUE
       CALL  IS02R (N, X, U, P, Y, Y2, DX, E, R, IERR)

Результаты:

       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