Текст подпрограммы и версий
iib9r_c.zip , iib9d_c.zip
Тексты тестовых примеров
tiib9r_c.zip , tiib9d_c.zip

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

Назначение

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

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

Пусть заданы узлы неравномepнoй прямоугольной сетки:

   X1A  =  { x1(1), x2(1), ..., xM(1) } ;   X2A  =  { x1(2), x2(2), ..., xN(2) } ,

 значения интерполируемой функции    f (x(1), x(2))

 и ее второй частной производной   ∂(2) f (x(1), x(2)) / ∂x(2) 2
 в узлах этой сетки:
 YA(I, J)  =  f (X1A(I), X2A(J)) ;   Y2A(I, J)  =  ∂2 f (X1A(I), X2A(J)) / ∂x(2) 2 ,

   I = 1, 2, ..., M ;  J = 1, 2, ..., N , 

а также следующие значения первой частной производной функции  f по первой переменной:

   YP1  =  ∂ f (X1A(1), X2)) / ∂x(1) ,   YPM  =  ∂ f (X1A(M), X2) / ∂x(1) . 

Здесь X2 - вторая координата заданной точки (X1, X2) внутри сеточного прямоугольника, в которой ищется значение интерполируемой функции: Y = f (X1, X2). Используется метод бикубических сплайнов.

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

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

    int iib9r_c (real *x1a, real *x2a, real *ya, real *y2a,
            integer *m, integer *n, real *x1, real *x2, real *yp1, real *ypm, 
            real *y, real *yt, real *y2t, real *yyt, real *u)

Параметры

x1a -
x2a  
вещественные векторы длины m и n, содержащие узлы  {x1 (1), x2 (1), ...,xm (1)} и {x1 (2), x2 (2), ...,xn (2)} соответственно;
ya - вещественный двумерный массив размеров m на n, содержащий значения интерполируемой функции  f (x (1), x (2)) двух переменных в узлах заданной сетки;
y2a - вещественный двумерный массив размеров m и n, содержащий значения второй частной производной функции  f по второй переменной в узлах заданной сетки;
m, n - длины векторов x1a и x2a соответственно (тип: целый);
x1, x2 - координаты заданной точки внутри сеточного прямоугольника, в которой ищется значение интерполируемой функции  f (тип: вещественный);
yp1 -
ypm  
заданные значения первой частной производной функции  f по первой переменной в точках (x1a (1), x2) и (x1a (m), x2) соответственно (тип: вещественный);
y - вещественная переменная, содержащая вычисленное значение интерполируемой функции  f в точке (x1, x2);
yt - вещественный вектор длины n, используемый в подпрограмме в качестве рабочего;
y2t - вещественный вектор длины  max (m, n), используемый в подпрограмме в качестве рабочего;
yyt, u - вещественные векторы длины m, используемые в подпрограмме в качестве рабочих.

Версии

i ib9d_c - выполняет те же вычисления, что и подпрограмма i ib9r_c, однако в режиме удвоенной точности. При этом все формальные параметры, кроме m и n, должны иметь тип double.

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

ids8r_c -
ids8d_c  
вычисление вторых производных таблично - заданной функции в узлах неравномерной сетки методом кубических или натуральных кубических сплайнов при заданных первых производных функции в концевых узлах сетки в режиме одинарной и удвоенной точности соответственно; используются в подпрограммах i ib9r_c и i ib9d_c соответственно;
i is8r_c -
i is8d_c  
вычисление значения в заданной точке интерполируемой табличной функции, определенной в узлах неравномерной сетки, упорядоченной по возрастанию, по известным значениям ее второй производной в узлах этой сетки методом кубических сплайнов в режиме одинарной и удвоенной точности соответственно; используются в подпрограммах i ib9r_c и i ib9d_c соответственно.

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

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

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

    /* Local variables */
    extern int iib9r_c(float *, float *, float *, float *, int *, int *,
                       float *, float *, float *, float *, float *, float *,
                       float *, float *, float *);
    static int i__, j, m, n;
    static float r__, u[5], y, x1, x2, ya[20] /* was [5][4] */,
               yt[4], x1a[5], x2a[4], y2a[20] /* was [5][4] */,
               yp1, y1t, y2t[5], ypm, yyt[5];
    int i__1, i__2;

#define ya_ref(a_1,a_2) ya[(a_2)*5 + a_1 - 6]
#define y2a_ref(a_1,a_2) y2a[(a_2)*5 + a_1 - 6]

    m = 5;
    n = 4;
    r__ = 0.f;
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        x1a[i__ - 1] = r__;
/* l1: */
        r__ += .1f;
    }
    r__ = 0.f;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        x2a[i__ - 1] = r__;
/* l2: */
        r__ += .1f;
    }
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            ya_ref(i__, j) = (float)sin(x1a[i__ - 1]) * (float)sin(x2a[j - 1]);
/* l4: */
            y2a_ref(i__, j) = -ya_ref(i__, j);
        }
    }
    x1 = .55f;
    x2 = .65f;
    yp1 = (float)cos(x1a[0]) * (float)sin(x2);
    ypm = (float)cos(x1a[m - 1]) * (float)sin(x2);
    iib9r_c(x1a, x2a, ya, y2a, &m, &n, &x1, &x2, &yp1, &ypm, &y, yt, y2t,
            yyt, u);
    y1t = (float)sin(x1) * (float)sin(x2);

    printf("\n %24.14e %24.14e \n",y, y1t);
    return 0;
} /* main */


Результат:   y  =  0.316370