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

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

Назначение

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

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

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

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

    int is02r_c (integer *n, real *x, real *u, real *p, real *y,
            real *y2, real *dx, real *e, real *r, integer *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_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы is02r_c. 
   Замечания по использованию
  Подпрограмма is02r_c может быть использована для сглаживания приближенно заданных на неравномерной сетке t1 < t2 < ...< tn значений  uk = u (tk),  vk = v (tk),  k = 1, 2, ..., n параметрической циклической функции  u = u (t),  v = v (t). Для этого достаточно дважды обратиться к подпрограмме is02r_c.

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

int main(void)
{
    /* Builtin functions */
    double sin(double);

    /* Local variables */
    extern int is02r_c(int *, float *, float *, float *, float *, float *,
                       float *, float *, float *, int *);
    static int ierr;
    static float e[2], h__;
    static int i__, n;
    static float p[9], r__[54] /* was [6][9] */, u[9], x[9], y[9], y2[9], pi,
                 dx[8], pi2;
    int i__1, i;

    pi = 3.141592f;
    pi2 = pi * 2;
    n = 9;
    h__ = pi2 / (n - 1);
    x[0] = 0.f;
    u[0] = 0.f;
    p[0] = 1.f;
    i__1 = n;
    for (i__ = 2; i__ <= i__1; ++i__) {
        x[i__ - 1] = x[i__ - 2] + h__;
        u[i__ - 1] = (float)sin(x[i__ - 1]);
        p[i__ - 1] = 1.f;
/* l1: */
    }
    is02r_c(&n, x, u, p, y, y2, dx, e, r__, &ierr);

    printf("\n %5i \n",ierr);
    for (i = 0; i <= 6; i += 3) {
        printf("\n %16.7e %16.7e %16.7e \n",y[i], y[i+1], y[i+2]);
    }
    for (i = 0; i <= 6; i += 3) {
        printf("\n %16.7e %16.7e %16.7e \n",y2[i], y2[i+1], y2[i+2]);
    }
    for (i = 0; i <= 4; i += 4) {
        printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                  dx[i], dx[i+1], dx[i+2], dx[i+3]);
    }
    printf("\n %16.7e %16.7e \n",e[0], e[1]);
    return 0;
} /* main */


Результаты:

       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