Текст подпрограммы и версий
idb8r_c.zip , idb8d_c.zip 		Тексты тестовых примеров
tidb8r_c.zip , tidb8d_c.zip

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

Назначение

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

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

Пусть заданы узлы неравномepнoй прямоугольной сетки:
X1A = { x1 ( 1 ), x2 ( 1 ),..., xM ( 1 )}; X2A = { x1 ( 2 ), x2 ( 2 ),..., xN ( 2 )}, значения функции  f (x ( 1 ), x ( 2 )) в узлах этой сетки: YA ( I, J ) = f (X1A ( I ), X2A( J )), I = 1, 2,..., M;  J = 1, 2,..., N, а также значения первой частной производной  ∂ f (x ( 1 ), x ( 2 )) / ∂ x ( 2 ) на нижней стороне сеточного прямоугольника:
YP1 ( I ) = ∂ f (X1A ( I ), X2A (1)) / ∂ x ( 2 ), I = 1, 2,..., M
и на верхней стороне сеточного прямоугольника:
YPN ( I ) = ∂ f (X1A ( I ), X2A (N)) / ∂ x ( 2 ),  I = 1, 2,..., M.

Подпрограмма   idb8r_c   вычисляет  вторые  частные  производные  ∂2 f (x ( 1 ), x ( 2 )) / ∂x ( 2 ) 2 в узлах заданной сетки методом одномерных кубических сплайнов по направлению  x ( 2 ).

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

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

    int idb8r_c (real *x1a, real *x2a, real *ya, integer *m,
            integer *n, real *yp1, real *ypn, real *y2a, real *yt, real *y2t, 
            real *rab)

   Параметры
x1a -
x2a   		вещественные векторы длины  m и  n, содержащие узлы { x1 ( 1 ), x2 ( 1 ),..., xm ( 1 )} и { x1 ( 2 ), x2 ( 2 ),..., xn ( 2 )} соответственно;
ya - вещественный двумерный массив размеров m на n, содержащий значения функции  f (x ( 1 ), x ( 2 )) двух переменных в узлах заданной сетки;
m, n - длины векторов x1a и x2a соответственно (тип: целый);
yp1 -
ypn  
 вещественные векторы длины m, содержащие значения первой частной производной функции  f по второй переменной на нижней и верхней сторонах сеточного прямоугольника соответственно;
y2a - вещественный двумерный массив размеров m на n, содержащий вычисленные значения второй частной производной функции  f по второй переменной в узлах заданной сетки;
yt, y2t -
rab   		вещественные векторы длины n, используемые в подпрограмме в качестве рабочих 
   Версии
idb8d_c - выполняет те же вычисления, что и idb8r_c, но в режиме удвоенной точности. 
   Вызываемые подпрограммы
       ids8r_c -
       ids8d_c   		вычисление вторых производных таблично - заданной функции в узлах неравномерной сетки методом кубических или натуральных кубических сплайнов при заданных первых производных функции в концевых узлах сетки в режимах одинарной и удвоенной точности соответственно; используются в подпрограммах idb8r_c и idb8d_c 
   Замечания по использованию:  нет

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

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

    /* Local variables */
    extern int idb8r_c(float *, float *, float *, int *, int *, float *,
                       float *, float *, float *, float *, float *);
    static int i__, j, m, n;
    static float r__, ya[20] /* was [5][4] */, yt[4], x1a[5], x2a[4],
                     y2a[20] /* was [5][4] */, yp1[5], y2t[4], rab[4],
                     ypn[5], y2at[20] /* was [5][4] */;
    int i__1, i__2;

#define ya_ref(a_1,a_2) ya[(a_2)*5 + a_1 - 6]
#define y2at_ref(a_1,a_2) y2at[(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__) {
        yp1[i__ - 1] = (float)sin(x1a[i__ - 1]) * (float)cos(x2a[0]);
/* l3: */
        ypn[i__ - 1] = (float)sin(x1a[i__ - 1]) * (float)cos(x2a[n - 1]);
    }
    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: */
            y2at_ref(i__, j) = -ya_ref(i__, j);
        }
    }
    idb8r_c(x1a, x2a, ya, &m, &n, yp1, ypn, y2a, yt, y2t, rab);

    for (i__ = 0; i__ <= 15; i__+= 5) {
        printf("\n %15.6e%15.6e%15.6e%15.6e%15.6e \n",
                 y2a[i__], y2a[i__+1], y2a[i__+2], y2a[i__+3], y2a[i__+4]);
    }
    for (i__ = 0; i__ <= 15; i__+= 5) {
        printf("\n %15.6e%15.6e%15.6e%15.6e%15.6e \n",
              y2at[i__], y2at[i__+1], y2at[i__+2], y2at[i__+3], y2at[i__+4]);
    }
    return 0;
} /* main */


Результаты:      

                     |   0.0                     0.0                    0.0                    0.0
                     |  -0.19949e-05   -0.99745e-02   -0.19551e-01   -0.29526e-01
       y2a  =    |  -0.37802e-05   -0.19849e-01   -0.39504e-01   -0.58755e-01
                     |  -0.57844e-05   -0.29525e-01   -0.58763e-01   -0.87398e-01
                     |  -0.85849e-05   -0.38906e-01   -0.77433e-01   -0.11517