Текст подпрограммы и версий
ip04r_c.zip , ip04d_c.zip
Тексты тестовых примеров
tip04r_c.zip , tip04d_c.zip

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

Назначение

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

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

Пусть X1L, X1U, X2L и X2U следующим образом определяют вершины элементарного сеточного прямоугольника:

        (X1L,X2U)   |                           |    (X1U,X2U)
      ______________________________________
                             |  4                   3  |
                             |       (X1,X2)       |

                             |  1                   2  |
      ______________________________________
        (X1L,X2L)   |                            |    (X1U,X2L)

Эти вершины пронумерованы на рисунке цифрами, начиная с нижней левой вершины против часовой стрелки. Пусть в каждой из вершин известны значения интерполируемой функции

        Y  =  f(x1, x2) ,
 обеих ее частных производных первого порядка:
        Y1  =  ∂f / ∂x1 ,   Y2  =  ∂f / ∂x2
 и смешанной производной   Y12  =  ∂2 f / ∂x1 ∂x2 . 

Подпрограмма ip04r_c при помощи бикубической интерполяции вычисляет в точке (X1,X2), заданной внутри исходного элементарного сеточного прямоугольника, значение функции

        AY = f (X1, X2) ,
 обеих ее частных производных первого порядка:
      AY1  =  ∂ f (X1, X2) / ∂x1 ,   AY2  =  ∂ f (X1, X2) / ∂x2
 и смешанной производной   AY12  =  ∂2 f (X1, X2) / ∂x1 ∂x2. 

Используются следующие формулы:

                          4     4                                               
           AY  =    ∑    ∑   ci j  t i-1 u j-1 ;
                         i=1  j=1
                          4     4
          AY1  =   ∑    ∑   (i - 1) ci j  t i-2 u j-1 ;
                         i=1  j=1
                          4     4
          AY2  =   ∑    ∑   (j - 1) ci j  t i-1 u j-2 ;
                         i=1  j=1
                          4     4
        AY12  =   ∑    ∑   (i - 1) (j - 1) ci j  t i-2 u j-2 ;
                         i=1  j=1 

Здесь  t = (X1 - X1L) / (X1U - X1L), U = (X2 - X2L) / (X2U - X2L). Коэффициенты  ci j  ( i, j = 1, 2, 3, 4 ) вычисляются в подпрограмме ip04r_c.

Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.

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

    int ip04r_c (real *y, real *y1, real *y2, real *y12, real *x1l,
             real *x1u, real *x2l, real *x2u, real *x1, real *x2, real *ay,
             real *ay1, real *ay2, real *ay12, real *c)

Параметры

          y, y1 -
        y2, y12  
вещественные векторы длины  4, содержащие соответственно значения интерполируемой функции, обеих ее частных производных первого порядка и смешанной производной в вершинах заданного элементарного сеточного прямоугольника, пронумерованных начиная с нижней левой вершины против часовой стрелки;
            x1l -
            x1u  
левая и правая координаты вершин элементарного сеточного пpямoyгольника по первой переменной (тип: вещественный);
            x2l -
            x2u  
нижняя и верхняя координаты вершин элементарного сеточного прямоугольника по второй переменной (тип: вещественный);
x1, x2 - координаты точки, заданной внутри исходного элементарного сеточного прямоугольника, в которой выполняется интерполяция (тип: вещественный);
      ay, ay1 -
       ay2, ay12  
вещественные переменные, значения которых полагаются равными соответственно значениям в точке (x1, x2) интерполируемой функции, ее частных производных первого порядка и смешанной производной;
c - вещественный двумерный массив размеров 4 на 4, в котором помещаются коэффициенты, используемые при бикубической интерполяции.

Версии

ip04d_c - бикубическая интерполяция функции двух переменных в режиме удвоенной точности; при этом все формальные параметры должны иметь тип double.

Вызываемые подпрограммы: нет

Замечания по использованию: нет

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

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

    /* Local variables */
    extern int ip04r_c(float *, float *, float *, float *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *);
    static float c__[16] /* was [4][4] */, y[4], x1, y1[4], y2[4], x2, ay,
            y12[4], ty, ay1, ay2, x1l, x2l, x1u, x2u, ty1, ty2, ay12, ty12;
    int i;

    x1l = 0.f;
    x1u = .1f;
    x2l = 0.f;
    x2u = .15f;
    y[0] = (float)sin(x1l) * (float)sin(x2l);
    y[1] = (float)sin(x1u) * (float)sin(x2l);
    y[2] = (float)sin(x1u) * (float)sin(x2u);
    y[3] = (float)sin(x1l) * (float)sin(x2u);
    y1[0] = (float)cos(x1l) * (float)sin(x2l);
    y1[1] = (float)cos(x1u) * (float)sin(x2l);
    y1[2] = (float)cos(x1u) * (float)sin(x2u);
    y1[3] = (float)cos(x1l) * (float)sin(x2u);
    y2[0] = (float)sin(x1l) * (float)cos(x2l);
    y2[1] = (float)sin(x1u) * (float)cos(x2l);
    y2[2] = (float)sin(x1u) * (float)cos(x2u);
    y2[3] = (float)sin(x1l) * (float)cos(x2u);
    y12[0] = (float)cos(x1l) * (float)cos(x2l);
    y12[1] = (float)cos(x1u) * (float)cos(x2l);
    y12[2] = (float)cos(x1u) * (float)cos(x2u);
    y12[3] = (float)cos(x1l) * (float)cos(x2u);
    x1 = .06f;
    x2 = .07f;
    ip04r_c(y, y1, y2, y12, &x1l, &x1u, &x2l, &x2u, &x1, &x2, &ay, &ay1, &ay2,
            &ay12, c__);

    for (i = 0; i <= 12; i+= 4) {
           printf("\n %16.7e %16.7e %16.7e %16.7e \n",
                   c__[i], c__[i+1], c__[i+2], c__[i+3]);
    }
    ty = (float)sin(x1) * (float)sin(x2);
    ty1 = (float)cos(x1) * (float)sin(x2);
    ty2 = (float)sin(x1) * (float)cos(x2);
    ty12 = (float)cos(x1) * (float)cos(x2);
           printf("\n %16.7e %16.7e \n",ay, ay1);
           printf("\n %16.7e %16.7e \n",ay2, ay12);
           printf("\n %16.7e %16.7e \n",ty, ty1);
           printf("\n %16.7e %16.7e \n",ty2, ty12);
    return 0;
} /* main */


Результаты:
        
           ay = 0.4194o5e-02 
         ay1 = 0.698169e-01 
         ay2 = 0.598171e-01 
       ay12 = 0.995755