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

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

Назначение

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

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

Вычисляются коэффициенты кубического сплайна 

 S (x) = yi + c1 i ( x - xi ) + c2 i ( x - xi )2 + c3 i ( x - xi )3 ,
 x ∈ [xi, xi + 1]  ,    i = 1, 2, ..., N - 1  ,

соответствующего сетке  {x i} и интерполирующего заданные значения табличной функции  yi = y (x i) ,  i = 1, 2, ..., N в узлах сетки.

Коэффициенты  c1 i , c2 i , c3 i ,  i = 1, 2, ..., N - 1 определяются значениями  S''i = S'' (x i) ,  i = 1, 2, ..., N, для нахождения которых решается система линейных алгебраических уравнений, полученная из условия непрерывности функции  S' (x) в узлах сетки и из краевых условий вида: 

          2 * S''1 + b1 * S''2 = b2
          b3 * S''N - 1 + 2 * S''N = b4  ,

где  b1 , b2 , b3 , b4 - заданные параметры. Если значения производных  y''1 = y'' (x1) и  y''N = y'' (xN) интерполируемой функции известны, то надо положить: 

          b1 = 0 ,  b2 = 2 * y''1
          b3 = 0 ,  b4 = 2 * y''N

Если известны значения  y'1 = y' (x1) и  y'N = y' (xN), то надо положить: 

 b1 = 1 ,   b2 = (6/(x2 - x1)) * ( (y2 - y1)/(x2 - x1) - y'1 )  ,
 b3 = 1 ,   b4 = (6/(xN - xN - 1)) * ( y'N - (yN - yN - 1)/(xN - xN - 1) )

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

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

    int iis3r_c (real *x, real *y, integer *nx, real *bpar,
            real *c, integer *ierr)

   Параметры
x - вещественный вектоp длины  nx, содержащий заданные значения узлов сетки и упорядоченный так, что  x (1) < x (2) < ... < x (nx);
y - вещественный вектоp длины  nx, содержащий заданные в узлах сетки значения интерполируемой функции;
nx - заданное число узлов сетки,  nx ≥ 2 (тип: целый);
bpar - вещественный вектоp длины 4, содержащий параметры краевых условий,  bpar (i) = bi,  i = 1, 2, 3, 4;
c - вещественный двумерный массив размера  3 * (nx - 1), содержащий вычисленные коэффициенты интерполяционного кубического сплайна;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
ierr=65 - когда заданное число узлов сетки меньше 2;
ierr=66 - когда элементы вектоpа  x не упорядочены по возрастанию. 
   Версии:  нет 
   Вызываемые подпрограммы
uti i10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы i is3r_с. 
   Замечания по использованию:  нет

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

int main(void)
{
    /* Initialized data */
    static float x[4] = { 0.f,.33f,1.f,2.f };

    /* Local variables */
    static float bpar[4];
    static int ierr;
    extern int iis3r_c(float *, float *, int *, float *, float *, int *);
    static float c__[9]  /* was [3][3] */;
    static int i__, j;
    static float y[4];
    static int nx;
    int i__1;

#define c___ref(a_1,a_2) c__[(a_2)*3 + a_1 - 4]

    nx = 4;
    i__1 = nx;
    for (i__ = 1; i__ <= i__1; ++i__) {
        y[i__ - 1] = ((x[i__ - 1] - 2.f) * x[i__ - 1] + 3.f) * x[i__ - 1] + 
                4.f;
/* l5: */
    }
    bpar[0] = 1.f;
    bpar[1] = 6.f / (x[1] - x[0]) * ((y[1] - y[0]) / (x[1] - x[0]) - 3.f);
    bpar[2] = 1.f;
    bpar[3] = 6.f / (x[3] - x[2]) * (7.f - (y[3] - y[2]) / (x[3] - x[2]));
    iis3r_c(x, y, &nx, bpar, c__, &ierr);

    printf("\n %5i \n",ierr);
/* l10: */
    for (j = 1; j <= 3; ++j) {
        printf("\n  %16.7e %16.7e %16.7e \n",
               c___ref(j-1, 0), c___ref(j-1, 1), c___ref(j-1, 2));
    }
    return 0;
} /* main */


Результаты:

      ierr = 0
 
                  |  3.0000  -2.0000  1.0000  |
      c__  =  |  2.0067  -1.0100  1.0000  |
                  |  2.0000   1.0000  1.0000  |