Участник:Без имени/Алгоритм параллельного вычисления интеграла с использованием составной формулы трапеций

Алгоритм параллельного вычисления интеграла с использованием составной формулы трапеций.

Основные авторы описания: Участник: Кузнецова Д. Д.

Свойства и структура алгоритма

Общее описание алгоритма

Метод трапеций является отдельным представителем класса алгоритмов численного интегрирования (приближённого вычисления определённых интегралов), также известных как квадратурные формулы. Алгоритм предназначен для решения задачи нахождения приближённого значения определённого интеграла от функции [math]\displaystyle{ f(x) }[/math] на отрезке [math]\displaystyle{ [a,b] }[/math]: [math]\displaystyle{ I = \int_a^b f(x)dx }[/math]

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

Требуется найти приближённое значение определённого интеграла от функции [math]\displaystyle{ f(x) }[/math] на отрезке [math]\displaystyle{ [a,b] }[/math]:

[math]\displaystyle{ I = \int_a^b f(x)dx }[/math]

Исходный отрезок [math]\displaystyle{ [a,b] }[/math] разбивается на [math]\displaystyle{ n }[/math] элементарных отрезков с помощью узлов интегрирования [math]\displaystyle{ x_i }[/math], где [math]\displaystyle{ i = 0,1,\dots,n }[/math], таких что:

[math]\displaystyle{ x_0 = a, x_n = b }[/math]

На каждом элементарном отрезке [math]\displaystyle{ [x_i,x_{i+1}] }[/math], где [math]\displaystyle{ i = 0,1,\dots,n-1 }[/math], применяется формула трапеций:

[math]\displaystyle{ \int_{x_i}^{x_{i+1}} f(x)dx \approx \frac{f(x_i) + f(x_{i+1})}{2} h_i }[/math]

Где [math]\displaystyle{ h_i = x_{i+1} - x_i }[/math] – длина элементарного отрезка. Суммируя приближённые значения интегралов по всем элементарным отрезкам, получаем составную формулу трапеций:

[math]\displaystyle{ \int_a^b f(x)dx \approx \sum_{i=0}^{n-1} \frac{f(x_i) + f(x_{i+1})}{2} h_i }[/math]

В частности, когда все элементарные отрезки имеют одинаковую длину ([math]\displaystyle{ h_i = h = \frac{b-a}{n} }[/math]), формула принимает вид:

[math]\displaystyle{ \int_a^b f(x)dx \approx I_h = \frac{h}{2} [f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)] }[/math]

[math]\displaystyle{ x_i = a + ih, i = 0,1,\dots,n }[/math]

Метод трапеций имеет второй порядок точности и является точным для многочленов степени не выше первой.

Вычислительное ядро алгоритма

Вычислительное ядро метода трапеций представляет собой вычисление суммы значений в узлах сетки:

[math]\displaystyle{ S = f(x_0) + f(x_1) + ... + f(x_n) }[/math]

Где

[math]\displaystyle{ f(x) }[/math] - подынтегральная функция

[math]\displaystyle{ x_i = a+ih }[/math] - узлы сетки

n - количество интервалов разбиения

[math]\displaystyle{ I \approx \frac{h}{2} [f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)] }[/math]

Основная вычислительная нагрузка приходится на [math]\displaystyle{ \sum_{i=1}^{n-1} f(x_i) }[/math]

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

Макроструктура алгоритма

1. Инициализация параметров

 INIT(a, b, n) → h, sum_initial

Вычисление шага сетки:

 h = (b - a) / n

Инициализация суммы:

 sum_initial = 0.5 * (f(a) + f(b))

2. Вычисление суммы значений функции

 SUM_FUNCTION_VALUES(a, h, n) -> sum_core

Основная вычислительная макрооперация

Вычисление: sum_core = ∑[i=1 to n-1] f(a + i·h)

3. Финальное вычисление интеграла

 FINALIZE(sum_initial, sum_core, h) -> integral

Объединение результатов:

 total_sum = sum_initial + sum_core

Финальное вычисление:

 integral = total_sum * h

Схема реализации последовательного алгоритма

1. Инициализация параметров интегрирования

// Параметры интегрирования
    double a = 0.0;    // нижний предел
    double b = 10.0;   // верхний предел
    int n = 1000000;   // количество интервалов

2. Вычисление шага сетки

double h = (b - a) / n;

- Вычисляется равномерное расстояние между соседними узлами

- Определяет плотность разбиения области интегрирования

3. Инициализация суммы граничных значений

double integral = 0.5 * (f(a) + f(b));

- Вычисляются значения функции на границах интервала

- Умножаются на коэффициент 0.5 согласно формуле трапеций

- Создается начальное приближение интеграла

4. Основной вычислительный цикл

for (int i = 1; i < n; i++) 
    {
        double x = a + i * h;
        integral += f(x);
    }

- Для каждого i от 1 до n-1 вычисляется x_i = a + i·h

- Для каждого узла вычисляется f(x_i). Это наиболее вычислительно затратная часть алгоритма

- Накопление суммы

5. Финальное вычисление интеграла

result = integral * h

Полная схема алгоритма

НАЧАЛО
  │
  |- Ввод параметров: a, b, n
  │
  |- Проверка: n > 0, a < b
  │
  |- Вычисление: h = (b - a) / n
  │
  |- Инициализация: sum = 0.5 * [f(a) + f(b)]
  │
  |- ЦИКЛ i = 1 ДО n-1:
  │   │
  │   |-Вычисление: x = a + i * h
  │   │
  │   |- Вычисление: fx = f(x)
  │   │
  │   |─ Накопление: sum = sum + fx
  │
   Результат: integral = sum * h

Последовательная сложность алгоритма

Метод трапеций относится к классу алгоритмов с линейной сложностью [math]\displaystyle{ O(n) }[/math] по количеству интервалов разбиения n. Основной вклад в сложность вносят вычисления функции [math]\displaystyle{ f(x) }[/math] именно от неё зависит количество проводимых вычислений в процессе применения формулы трапеций от каждого участка.

Информационный граф

V1 - получение параметров интегрирования, V2.i - рассылка границ на каждый процесс, V3.i - вычисление площадей на каждом отрезке, V4 - итоговое суммирование

Рис. 1. - Информационный граф, V1 - получение параметров интегрирования, V2.i - рассылка границ на каждый процесс, V3.i - вычисление площадей на каждом отрезке, V4 - итоговое суммирование

Ресурс параллелизма алгоритма

Основная вычислительная нагрузка сосредоточена в полностью независимых операциях вычисления значений функции в узлах сетки. Параллельная сложность составляет [math]\displaystyle{ O(1) }[/math] для основной расчетной части при использовании достаточного количества процессоров. Структура параллелизма обеспечивает балансировку нагрузки и простоту реализации. При больших n доля последовательных операций стремится к нулю.

Входные и выходные данные алгоритма

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

Программная реализация алгоритма

Особенности реализации последовательного алгоритма

Для большинства практических применений метода трапеций достаточно одинарной точности (float). В большинстве практических случаев не требует двойной точности вычислений, поскольку его основная погрешность определяется теоретической ошибкой метода O(1/n²), а не вычислительной ошибкой округления.

Возможные способы и особенности параллельной реализации алгоритма

Основной способ распределения данных. Весь диапазон [a, b] разбивается на p блоков (где p - количество процессов), каждый процесс получает свой непрерывный блок интервалов. Для балансировки нагрузки при неравномерном количестве интервалов применяется циклическое распределение остатка - первые r процессов получают на один интервал больше.

Ключевая коммуникационная операция - глобальная редукция (MPI_Reduce) для суммирования частичных интегралов. Объем пересылаемых данных минимален - каждый процесс отправляет всего одно вещественное число (свой частичный интеграл).

Параллельная эффективность ограничена последовательной долей по закону Амдала, которая включает инициализацию и финальную редукцию. При больших n эта доля становится пренебрежимо малой. Для достижения максимальной производительности рекомендуется выбирать размер задачи n таким, чтобы каждый процесс/поток получал достаточно большую порцию работы (обычно тысячи интервалов), минимизируя относительные накладные расходы.

Результаты прогонов

Processes: 2, Time: 0.080949 seconds

Processes: 2, Time: 0.074890 seconds

Processes: 2, Time: 0.058663 seconds

Processes: 2, Time: 0.058114 seconds

Processes: 4, Time: 0.047105 seconds

Processes: 4, Time: 0.043581 seconds

Processes: 4, Time: 0.048686 seconds

Processes: 4, Time: 0.046107 seconds

Processes: 8, Time: 0.024348 seconds

Processes: 8, Time: 0.022984 seconds

Processes: 8, Time: 0.021923 seconds

Processes: 8, Time: 0.024774 seconds

Processes: 16, Time: 0.013365 seconds

Processes: 16, Time: 0.011646 seconds

Processes: 16, Time: 0.011226 seconds

Processes: 16, Time: 0.010894 seconds

Processes: 20, Time: 0.008676 seconds

Processes: 20, Time: 0.009282 seconds

Processes: 20, Time: 0.009332 seconds

Processes: 20, Time: 0.010975 seconds

Processes: 24, Time: 0.008456 seconds

Processes: 24, Time: 0.008577 seconds

Processes: 24, Time: 0.008664 seconds

Processes: 24, Time: 0.009281 seconds

Processes: 28, Time: 0.007415 seconds

Processes: 28, Time: 0.007020 seconds

Processes: 28, Time: 0.007767 seconds

Processes: 28, Time: 0.006456 seconds

Литература

[1] Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. – 432 с.

[2] Якобовский М.В. Введение в параллельные методы решения задач, 2013. - 328 с.

[3] Антонов А.С. Технологии параллельного программирования MPI и OpenMP, 2012. - 344 с.