Текст подпрограммы и версий
qts8r_c.zip  qts8d_c.zip  		Тексты тестовых примеров
tqts8r_c.zip  tqts8d_c.zip 

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

Назначение

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

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

Вычисляется значение двойнoгo oпpеделенного интеграла от табличной функции f (x, y), заданной на неравномерной сетке (xi, yj),   i = 1, ..., N,  j = 1, ..., M, по квадратурной формуле, точной для многочленов второй степени.

Н.С.Бахвалов. Численные методы, "Hаука", M., 1975.

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

    int qts8r_c (real *rint, real *x1, real *x2, integer *n,
            integer *m, real *f, real *b1)

   Параметры
rint - вещественная переменная, содержащая вычисленное значение интеграла;
x1 - вещественный вектоp длины n, содержащий узлы неравномерной сетки по x;
x2 - вещественный вектоp длины m, содержащий узлы неравномерной сетки по y;
n - заданное число узлов сетки по x (тип: целый);
m - заданное число узлов сетки по y (тип: целый);
f - вещественный двумерный массив размера n на m, содержащий значения функции f (x, y);
b1 - вещественный вектоp длины n, используемый в подпрограмме в качестве рабочего. 
   Версии
qts8d_c - вычисление с удвоенной точностью определенного двухкратного интеграла от табличной функции, заданной на неравномерной сетке, по формуле типа Симпсона. 

   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

Требуется, чтобы  xi > xi - 1,  i = 2, ..., n,  yj > yj - 1,  j = 2, ..., m;  n ≥ 3,  m ≥ 3.

В подпрограмме qts8d_c параметры rint, x1, x2, f, b1 имеют тип double.

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

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

    /* Local variables */
    static float rint, f[2000]   /* was [50][40] */;
    static int i__, j, m, n;
    extern int qts8r_c(float *, float *, float *, int *,
                       int *, float *, float *);
    static float x, b1[50], h1, h2, x1[50], x2[40];
    int i__1, i__2;

#define f_ref(a_1,a_2) f[(a_2)*50 + a_1 - 51]

    h1 = .064114081632653058f;
    h2 = .08055358974358974f;
    n = 50;
    m = 40;
    x1[0] = 0.f;
    x2[0] = 0.f;
    x1[1] = h1 / 4.f;
    x2[1] = h2 / 4.f;
    i__1 = n;
    for (i__ = 3; i__ <= i__1; ++i__) {
/* l1: */
        x1[i__ - 1] = x1[i__ - 3] + h1 * 2.f;
    }
    i__1 = m;
    for (j = 3; j <= i__1; ++j) {
/* l2: */
        x2[j - 1] = x2[j - 3] + h2 * 2.f;
    }
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = m;
        for (j = 1; j <= i__2; ++j) {
            x = x1[i__ - 1] + x2[j - 1];
/* l3: */
            f_ref(i__, j) = (float)sin(x);
        }
    }
    qts8r_c(&rint, x1, x2, &n, &m, f, b1);

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


Результат:

       rint  =  0.21648