Текст подпрограммы и версий
qtshr_c.zip  qtshd_c.zip  		Тексты тестовых примеров
tqtshr_c.zip  tqtshd_c.zip 

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

Назначение

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

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

Вычисляются значения неопрeдeлeнного двойного интеграла по прямоугольникам x1 ≤ x ≤ xi,  y1 ≤ y ≤ yj,   i = 1, 2, ..., N,  j = 1, 2, ..., M от табличной функции f (x, y), заданной на равномерной сетке xi = x1 + (i - 1) h1,   yj = y1 + (j - 1) h2, по формуле Симпсона.

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

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

    int qtshr_c (real *rint, real *h1, real *h2, real *f,
            integer *n, integer *m)

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

Требуется, чтобы  n ≥ 4  и  m ≥ 4.

В подпрограмме qtshd_c параметры rint, h1, h2, f имеют тип double.

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

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

    /* Local variables */
    static float rint[2000] /* was [50][40] */,
                    f[2000] /* was [50][40] */;
    static int i__, j, m, n;
    static float x, y, h1, h2;
    extern int qtshr_c(float *, float *, float *, float *, int *, int *);

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

    h1 = .064114081632653058f;
    h2 = .08055358974358974f;
    n = 50;
    m = 40;
    y = 0.f;
    for (j = 1; j <= 40; ++j) {
        x = y;
        for (i__ = 1; i__ <= 50; ++i__) {
            f_ref(i__, j) = (float)sin(x);
/* l1: */
            x += h1;
        }
/* l2: */
        y += h2;
    }
    qtshr_c(rint, &h1, &h2, f, &n, &m);

        printf("\n  %12.6f %12.6f %12.6f %12.6f \n",
            rint_ref(1, 1), rint_ref(2, 2), rint_ref(n-1, m), rint_ref(n, m));
    return 0;
} /* main */


Результаты:

       rint(1, 1)         =  0.00000
       rint(2, 2)         =  0.00037
       rint (n-1, m)  =  0.12815
       rint (n, m)      =  0.00001