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

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

Назначение

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

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

Пусть на равномерной прямoyгoльной сетке (xi, yj),  i = 1, ..., M1,  j = 1, ..., M2, с шагами H1 и H2 по  x и  y соответственно, заданы значения функции  fi j = f (xi, yj).

По точке (X Y), в которой требуется произвести вычисления, определяются номеpа K и L такие, что X ∈ [xK, xK + 1],  Y ∈ [yL, yL + 1]. Значение производной порядка P по  x и  порядка Q по  y  ( 0 ≤ P ≤ 2N - 2, 0 ≤ Q ≤ 2N - 2 ) при K, L:  N ≤ K ≤ M1 - N, N ≤ L ≤ M2 - N вычисляется по формуле: 

      S (P,Q) (X,Y)  =  CN2 * 
                                               K+N         L+N       
                                         *     ∑            ∑           fi j s(P)i,2N (X)  s(Q)j,2N (Y) ,
                                             i=K-N+1    j=L-N+1

где C1 = 1,  C2 = 2/3,  C3 = 11/20,  C4 = 151/315,  si, 2N (Х) и  sj, 2N (Y) - суть  В - сплайны порядка 2N, а N = 1, 2, 3, 4. При других K и L вычисления ведутся по формуле Тейлора.

Гребенников А.И. Алгоритм быстрой аппроксимации функций двух переменных и их производных сплайнами на равномерной сетке. B сб. "Численный анализ на ФОРТРАНе", вып.20, Изд-во МГУ, 1977.

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

    int id24r_c (integer *p, integer *q, integer *n1, integer *m1,
            integer *m2, real *h1, real *h2, real *a1, real *a2, real *x, real *y,
             real *f, real *s, integer *ierr)

   Параметры
p, q - целые переменные, значения которых на единицу больше порядка вычисляемой производной по  x и  y соответственно (при p = 1, q = 1 вычисляется значение сплайна) 1 ≤ p ≤ 2n - 1, 1 ≤ q ≤ 2n - 1;
n1 - порядок сплайна, n1 = 2n,  1 ≤ n ≤ 4 (тип: целый);
m1 - число узлов на оси  x (тип: целый);
m2 - число узлов на оси  y (тип: целый);
h1 - шаг сетки по оси  x (тип: вещественный);
h2 - шаг сетки по оси  y (тип: вещественный);
a1 - координата первого узла на оси  x (тип: вещественный);
a2 - координата первого узла на оси  y (тип: вещественный);
x, y - координаты точки, в которой производятся вычисления производной (тип: вещественный);
f - двумерный вещественный массив размерности m1 на m2 заданных значений функций;
s - вещественная переменная, содержащая результат вычислений;
ierr - целая переменная, значение которой указывает причину возможной ошибки:
ierr=1 - шаги сетки слишком малы;
ierr=2 - точка (x, y) лежит вне допустимой области [x1, xm1] * [y1, ym2];
ierr=3 - невеpно задан порядок сплайна;
ierr=4 - невеpно задан порядок производной. 
   Версии:  нет 
   Вызываемые подпрограммы
uti i12_c - подпрограмма выдачи диагностических сообщений. 
   Замечания по использованию:  нет

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

int main(void)
{
    /* Builtin functions */
    double exp(double);

    /* Local variables */
    extern int id24r_c(int *, int *, int *, int *, int *, float *, float *,
                       float *, float *, float *, float *, float *, float *,
                       int *);
    static int ierr;
    static float f[121]  /* was [11][11] */;
    static int i__, j;
    static float s, x[11], y[11];
    static int ip, iq;

#define f_ref(a_1,a_2) f[(a_2)*11 + a_1 - 12]

    ip = 1;
    iq = 1;
    for (i__ = 1; i__ <= 11; ++i__) {
        x[i__ - 1] = (i__ - 1) * .1f;
        for (j = 1; j <= 11; ++j) {
            y[j - 1] = (j - 1) * .1f;
/* l5: */
            f_ref(i__, j) = (float)exp((float)(-x[i__ - 1] * x[i__ - 1] -
                    y[j - 1] * y[j - 1]));
        }
    }
    id24r_c(&ip, &iq, &c__4, &c__11, &c__11, &c_b5, &c_b5, &c_b7, &c_b7, &c_b9,
             &c_b10, f, &s, &ierr);

    printf("\n %16.7e \n",s);
    return 0;
} /* main */


Результатат:   s = 0.70865