1. Полученные алгоритмы легко обобщаются на случай многомерных пространств. Это осуществляется путем генерации в функции чисел, отвечающих за нужную координату. Так, в приложениях 3, 4 приведены тексты функций, реализующих вычисление двойных интегралов на последовательной и параллельной ЭВМ.
2. Для вычисления n-кратного интеграла необходимо знать величину n-1-кратного интеграла. Поэтому при расчете интегралов высокой кратности необходимо рассчитывать распределение вычислений по процессам особенно тщательно. Для решения этой проблемы в MPI предусмотрены виртуальная топология, группы, коммуникаторы.
3. Вследствие низкой сходимости метода Монте-Карло целесообразно применять его только при расчетах интегралов кратности 2 и выше.
ЛИТЕРАТУРА
1. Бахвалов Н.С. Численные методы. – М.: «Наука», 1975. – 631 с.
2. Грэхем Р., Кнут Д., Паташник О. Конкретная математика. – М.: «Мир», 1998. – 703 с.
3. Демидович Б.П. Сборник задач и упражнений по математическому анализу. – М.: «Наука», 1972. – 544 с.
4. Демидович Б.П., Марон И.А. Основы вычислительной математики. – М.: «Наука», 1966. – 664 с.
5. Информатика. Энциклопедический словарь для начинающих. – М.: «Педагогика-пресс», 1994. – 349 с.
6. Костевич Л.С., Лапко А.А. Теория игр. Исследование операций. – Мн.: «Вышэйшая школа», 1982. – 210 с.
7. Куст Ю., Юмагулов М. Математика. Основы математического анализа. – М.: «Айрис Пресс», 1999. – 270 с.
8. Шпаковский Г.И. Организация параллельных ЭВМ и суперскалярных процессоров. – Мн.: БГУ, 1996. – 283 с.
9. Шпаковский Г.И., Серикова Н.В. Программирование для многопроцессорных систем в стандарте MPI. – Мн.: БГУ, 2002. – 324 с.
10. Что такое Beowulf? (http://acdwp.narod.ru/la/podr/05020010.htm).
ПРИЛОЖЕНИЕ 1
// Файл DataType.h
// Содержит определенные и используемые в работе типы данных
// 20.05.04 г.
// Целые знаковые и беззнаковые
typedef unsigned int CODE_ERR; // Тип для возврата результата-ошибки
typedef unsigned int UINT; // Тип для счетчиков
typedef unsigned long int ULONG; // Тип для счетчиков
// Вещественные знаковые и беззнаковые
typedef float PTCOUNT; // Хранит кол-во точек
typedef long double LIMIT; // Тип для границ диапазона
typedef long double FUNC_TYPE; /* Тип результата подынтегральной функции */
typedef long double ARG_TYPE; // Тип для аргументов функций
typedef unsigned long double UFUNC_TYPE; // Тип для площади
typedef unsigned long double ERR_T; // Тип для точности
// Указатели на функции
typedef FUNC_TYPE (*SURF) (ARG_TYPE, ARG_TYPE); /* Функция двух аргументов */
typedef FUNC_TYPE (*CURVLINE) (ARG_TYPE); // Функция одного аргумента.
ПРИЛОЖЕНИЕ 2
// __MyHeader__.h
// Инлайн реализации вспомогательных функций
// Последние изменения 20 мая 2004 г.
///////////////////////////////////////////////////////////////////////////////////////////////////////////
// Включения
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include "DataType.h"
///////////////////////////////////////////////////////////////////////////////////////////////
inline CODE_ERR TotTimePrint
(clock_t start, clock_t stop, FILE *fp, char mode = 'f')
/* Инлайн функция печатает в файл (по умолчанию) или на дисплей (mode == 'd') время в минутах и секундах между отметками start и stop, полученными с помощью функции clock(). При этом, если время меньше 1 секунды функция выводит, что общее время выполнения чего-либо – < 1 sec. Возвращает код ошибки: "0" – если ошибок нет, "−1" – если значения start и/или stop – ошибочные или указан несуществующий режим */
{
if (start < 0 || stop < 0) return −1;
// Подсчитываем время выполнения
clock_t t = (stop-start) / CLK_TCK;
switch (mode) {
case 'f': // Печать в файл
if (t < 1) fprintf (fp, "Total time < 1 sec\n");
else fprintf (fp, "Total time = %d sec (%d min %d
sec)\n", t, t/60, (t - (t/60)*60));
break;
case 'd': // Печать на экран
if (t < 1) printf ("Total time < 1 sec\n");
else printf ("Total time = %d sec (%d min %d sec)\n",
t, t/60, (t - (t/60)*60));
break;
default:
return −1;
break;
}
return 0;
}
inline long double random (LIMIT left, LIMIT right)
/*Инлайн функция возвращает случайное число в интервале от left до right */
{ return ((long double)rand()) * (right-left) / RAND_MAX + left; }
ПРИЛОЖЕНИЕ 3
// Файл IntegralMK.h
/* Описания функций для подсчета интегралов методом Монте-Карло */
// 4 мая 2004 г. – 25 мая 2004 г.
///////////////////////////////////////////////////////////////////////////////////////////////
// Включения
#include <mpi.h>
#include "DataType.h"
#include "__MyHeader__.h"
///////////////////////////////////////////////////////////////////////////////////////////////
// Константы
#define EPS_ERR 1 // Неккоректная величина абсолютной ошибки
#define TOTPT_MAX 1e30 // Максимально допустимое число точек
///////////////////////////////////////////////////////////////////////////////////////////////
// Поверхность
FUNC_TYPE Surf (ARG_TYPE x, ARG_TYPE y);
// Границы у-правильной области, т.е. bot(x) и top(x)
FUNC_TYPE bot (ARG_TYPE arg); // Нижняя граница
FUNC_TYPE top (ARG_TYPE arg); // Верхняя граница
// Подынтегральная функция для 1-кратного интеграла
FUNC_TYPE f (ARG_TYPE arg);
///////////////////////////////////////////////////////////////////////////////////////////////
// Интерфейсы функций
FUNC_TYPE Integral1DSF
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-4);
// Рассчет 1-кратного интеграла по формуле парабол (Симпсона)
/* Использовалась для сравнения со значениями, полученными методом Монте-Карло */
// П Р О Т Е С Т И Р О В А Н А ! ! ! – 22.05.04 г.
FUNC_TYPE Integral1DMK
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3,
PTCOUNT totpt = 100, ULONG mult = 2);
// Вычисление 1-кратного интеграла методом Монте-Карло
// Вариант со средним значением функции
// П Р О Т Е С Т И Р О В А Н А ! ! ! – 25.05.04 г.
FUNC_TYPE Integral1DMK_par
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3,
PTCOUNT totpt = 100, ULONG mult = 2);
// Вычисление 1-кратного интеграла методом Монте-Карло на кластере MPI
// Вариант со средним значением функции
// П Р О Т Е С Т И Р О В А Н А
// НА КЛАСТЕРЕ ! ! ! – 26.05.04 г.
FUNC_TYPE Integral1DMK_par1
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3,
PTCOUNT totpt = 100, ULONG mult = 2);
/* Вычисление 1-кратного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции. Каждый процесс работает в более узком диапазоне, чем первоначальный */
// П Р О Т Е С Т И Р О В А Н А
// НА КЛАСТЕРЕ ! ! ! – 26.05.04 г.
FUNC_TYPE Integral1DMK_par2
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-3,
PTCOUNT totpt = 100, ULONG mult = 2);
/* Вычисление 1-кратного интеграла методом Монте-Карло на кластере MPI. Вариант со средним значением функции. Каждый процесс работает в более узком диапазоне, чем начальный. Операция редукции проводится однажды, когда каждый процесс подсчитал интеграл на своем узком участке */
// П Р О Т Е С Т И Р О В А Н А
// НА КЛАСТЕРЕ ! ! ! – 26.05.04 г.
UFUNC_TYPE Square
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps = 1e-4);
/* Вычисление площади под кривой f(x). Используется формула Симпсона */
// П Р О Т Е С Т И Р О В А Н А ! ! ! – 24.05.04 г.
inline UFUNC_TYPE RegionY
(CURVLINE bottom, CURVLINE top, LIMIT left,
LIMIT right, unsigned char var, ERR_T eps = 1e-4)
/* Вычисление площади между кривыми bottom(x) и top(x). Кривые сшиваются в точках left и right */
// П Р О Т Е С Т И Р О В А Н А ! ! ! – 25.05.04 г.
{
// var определяет случай вычисления площади области
switch (var) {
// Кривые top и bottom лежат в положительной полуоси
case 1: return Square(top, left, right, eps) –
Square(bottom, left, right, eps);
break;
// Кривые top и bottom лежат в отрицательной полуоси
case 2: return Square(bottom, left, right, eps) –
Square(top, left, right, eps);
break;
// Кривая top лежит в положительной полуоси
// кривая bottom лежит в отрицательной полуоси
case 3: return Square(top, left, right, eps) +
Square(bottom, left, right, eps);
break;
default: return 0;
}
}
FUNC_TYPE Integral2DMK
(SURF surf, LIMIT left, LIMIT right,
CURVLINE bottom, CURVLINE top, unsigned char var,
ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2);
/* Рассчет двойного интеграла для у-правильной области. Вариант со средним значением функции */
FUNC_TYPE Integral2DMK_par
(SURF surf, LIMIT left, LIMIT right,
CURVLINE bottom, CURVLINE top, unsigned char var,
ERR_T eps = 1e-3, PTCOUNT totpt = 100, ULONG mult = 2);
/* Рассчет двойного интеграла для у-правильной области на кластере MPI. Вариант со средним значением функции */
ПРИЛОЖЕНИЕ 4
// Файл IntegralMK.h
// Реализация метода Монте-Карло для подсчета интегралов
// Начало – 4 мая 2004 г.
///////////////////////////////////////////////////////////////////////////////////////////////
#include <math.h>
#include <limits.h>
#include "IntegralMK.h"
///////////////////////////////////////////////////////////////////////////////////////////////
FUNC_TYPE Integral1DSF
(CURVLINE f, LIMIT left, LIMIT right, ERR_T eps)
// Реализация формулы Симпсона для однократных интегралов
/* Использовалась для сравнения со значениями, полученными методом Монте-Карло */
// П Р О Т Е С Т И Р О В А Н А ! ! ! – 22.05.04 г.
{
if (eps <= 0) exit (EPS_ERR);
FUNC_TYPE _s, s = 0.0;
FUNC_TYPE t0 = f(left) + f(right);
ULONG n = 1024;
do {
_s = s;
ARG_TYPE h = (right - left) / (2*n);
FUNC_TYPE t1 = 0, t2 = 0;
// t1 - сумма нечетных, t2 – сумма четных
for (ULONG j = 1; j <= 2*n-1; j++)
(j % 2 == 0) ? t2 += f(left+j*h) : t1 += f(left+j*h);
s = (t0 + 4*t1 + 2*t2) *(h/3);
n *= 2; // mult == 2
} while (fabsl(s – _s) > eps && n <= ULONG_MAX);
return s;
}
///////////////////////////////////////////////////////////////////////////////////////////////
FUNC_TYPE Integral1DMK (CURVLINE f, LIMIT left, LIMIT right,
ERR_T eps, PTCOUNT totpt, ULONG mult)
// Вычисление однократного интеграла методом Монте-Карло.
// Вариант со средним значением функции
// П Р О Т Е С Т И Р О В А Н А ! ! ! – 25.05.04 г.
{
// Проверка корректности введенных данных
if (eps <= 0 || mult < 1 || totpt < 10 || totpt >= TOTPT_MAX) exit(1);
FUNC_TYPE _s = 0.0, s;
ARG_TYPE h, left_lim;
clock_t start = clock(), stop; // Запускаем таймер
do {
s = _s; // Запоминаем предыдущее значение интеграла
_s = 0.0; // Обнуляем сумму
// Инициализируем генератор случайных чисел временем
srand((UINT)(time(NULL)));
do { // Обеспечиваем малый шаг
h = (right - left) / totpt;
if (h > 1) totpt *= mult;
} while (h > 1);
if (!h) exit(1);
left_lim = left;
// Накапливаем сумму
while ((h > 0 && left_lim + h <= right) ||
(h < 0 && left_lim + h >= right)) {
_s += f(random(left_lim, left_lim+h));
left_lim += h;
}
// Вычисляем интеграл по формуле (3.1.10)
_s = (_s/totpt)*(right-left);
// Вывод на экран промежуточных данных