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

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

Назначение

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

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

Hа заданном элементарном пpямoyгoльнике  R = [xlx, xlx + 1] * [yly, yly + 1] вычисляются коэффициенты  CRm n натурального бикубического сплайна 

                          4       4            
        S (x, y)  =  ∑     ∑   CRm n (x - xlx)m -1 (y - yly)n - 1 ,
                       m = 1  n = 1
               (x, y) ∈ R  ,   1 ≤ lx ≤ NX  ,   1 ≤ ly ≤ NY  ,

соответствующего двумерной сетке {x i, yj} и интерполирующего заданные значения  fi j = f (x i, yj),  i = 1, 2, ..., NX,   j = 1, 2, ..., NY, в узлах сетки. Для нахождения 16 - ти коэффициентов  CRm n,  m, n = 1, 2, 3, 4, решается система линейных алгебраических уравнений, полученная по значениям бикубического сплайна в 16 - ти точках заданного прямоугольника.

Эти значения вычисляет библиотечная подпрограмма i is5r_c.

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

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

    int iis6r_c(real *f, real *x, integer *nx, real *y,
        integer *ny, integer *lx, integer *ly, real *c, real *wk,
        integer *ierr)

   Параметры
f - вещественный двумерный массив размера  nx * ny, содержащий заданные в узлах двумерной сетки значения интерполируемой функции;
x - вещественный вектоp длины  nx, содержащий заданные значения узлов одномерной сетки по оси  x и упорядоченный так, что  x (1) < x (2) < ... < x (nx);
nx - заданное число узлов одномерной сетки по оси  x;   nx ≥ 2 (тип: целый);
y - вещественный вектоp длины  ny, содержащий заданные значения узлов одномерной сетки по оси  y и упорядоченный так, что  y (1) < y (2) < ... < y (ny);
ny - заданное число узлов одномерной сетки по оси  y;   ny ≥ 2 (тип: целый);
lx - заданное число, определяющее элементарный прямоугольник  R по оси  x,  1 ≤ lx < nx (тип: целый);
ly - заданное число, определяющее элементарный прямоугольник  R по оси  y,  1 ≤ ly < ny (тип: целый);
c - вещественный двумерный массив размера 4 * 4, содержащий вычисленные на заданном элементарном прямоугольнике  R = [x lx, x lx + 1] * [yly, yly + 1] коэффициенты интерполяционного бикубического сплайна;
wk - вещественный вектоp длины  4 * max (4, ny) + max ((nx - 1) * 3, (ny - 1) * 3 + 4), используемый подпрограммой как рабочий;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
ierr=65 - когда заданное число узлов сетки по  x меньше 2;
ierr=66 - когда заданное число узлов сетки по  y меньше 2;
ierr=67 - когда элементы вектоpа  x не упорядочены по возрастанию;
ierr=68 - когда элементы вектоpа  y не упорядочены по возрастанию;
ierr=69 - когда  lx ≥ nx или  lx < 1;
ierr=70 - когда  ly ≥ ny или  ly < 1. 
   Версии:  нет 
   Вызываемые подпрограммы
i is5r_c - вычисление значений интерполяционного натурального бикубического сплайна на заданном множестве точек.
i is4r_c - вычисление значений интерполяционного кубического сплайна в заданных точках.
i is3r_c - интерполяция кубическими сплайнами табличной функции от одной переменной на неравномерной сетке при известных краевых условиях на вторые производные.
uti i10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы i is6r_c. 
   Замечания по использованию:  нет

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

int main(void)
{
    /* Initialized data */
    static float x[3] = { 1.f,3.f,5.f };
    static float y[3] = { 1.f,3.f,5.f };
    static float f[9] /* was [3][3] */ = { 2.f,4.f,3.f,.5f,1.f,.75f,1.f,2.f,
                                          1.5f };

    /* Local variables */
    static int ierr;
    extern int iis6r_c(float *, float *, int *, float *, int *, int *,
                       int *, float *, float *, int *);
    static float c__[16] /* was [4][4] */;
    static int i__;
    static float wk[26];
    static int lx, ly, nx, ny;

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

    nx = 3;
    ny = 3;
    lx = 1;
    ly = 1;
    iis6r_c(f, x, &nx, y, &ny, &lx, &ly, c__, wk, &ierr);

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


Результаты:

       ierr  =  0

                |  2.000000  -1.000000   0.000000    0.062500  |
    c__  =  |  1.375000  -0.687500   0.000000    0.042969  |
                |  0.000000   0.000000   0.000000    0.000000  |
                | -0.093750   0.046875   0.000000  -0.002930  |

   т.е. для  (x, y) ∈ [1., 3.] * [1., 3.]

                          4       4            
        S (x, y)  =  ∑     ∑   Cm n (x - 1.)m -1 (y - 1.)n - 1 
                       m = 1  n = 1