Участник:Greno4ka/Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки

Авторы: Устинов Г.А.(п. 1.6-1.10, 2.4, 2.7) Лиджиев Д.В.(п. 1.1-1.6)

1 ЧАСТЬ. Свойства и структура алгоритмов

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

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

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

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

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

Пусть на отрезке [math] [a,b] [/math] задана равномерная сетка, содержащая [math] n+1 [/math] узел:

[math] x_{i} = a + \frac{(b - a)} {n} i, i = 0,...,n [/math]

Тогда, согласно методу трапеций, можно численно найти определенный интеграл от функции [math] f (x) [/math] на отрезке [math][a,b][/math]:

[math] \int_a^b{f(x)dx} = \frac{b - a}{n} \left( \frac{f(x_0)}{2} + \sum^{n-1}_{i=1} {f(x_i)} + \frac{f(x_n)}{2} \right)[/math]

Будем полагать значение I найденным с точностью ε, если выполнено условие

[math] |I_{n_1} - I_{n_2}| \leq \epsilon |I_{n_2}|, n_1 \gt n_2, [/math] где [math] n_1 [/math] и [math] n_2 [/math] - количество узлов, образующих две разные сетки.

Данный метод неэффективен в силу ряда причин:

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

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

Вычислительное ядро алгоритма состоит из основного цикла алгоритма и сложность зависит от количества проходов по нему. То есть основное время работы алгоритма приходится на следующие операции:
- нахождение середины отрезка и значения функции в этой точке;
- вычисление частичных сумм на половинах отрезка;
- работа со стеком.

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

Алгоритм итерационный. Сначала значение текущего отрезка равно общему отрезку интегрирования. На каждой итерации выполняются следующие действия:
- На текущем отрезке считается приближенное значение интеграла
- Текущий отрезок делится пополам, и на получившихся отрезках высчитываются интегралы
- Сравнивается интеграл на всём отрезке с суммой двух интегралов на половинах. Если достигнута заданная точность, то:

1. Добавить значение интеграла текущего отрезка к общему интегралу
   
2. Если стек пуст, то завершить алгоритм
   
3. Иначе достать из стека отрезок и присвоить его значение текущему отрезку.
   

- Иначе

1. Одну половину отрезка положить в стек
   
2. Текущему отрезку присвоить значение второй половины
   

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

struct SumRecord {
	double A;
	double B;
	double s;
};

double JLS(double A, double B) {
	stack<SumRecord> localStack;
	double J = 0;
	double sAB = (f(A) + f(B)) * (B - A) / 2;
	SumRecord tmpRecordAB;
	while (1) {
		double C = (A + B) / 2;
		double sAC = (f(A) + f(C)) * (C - A) / 2;
		double sCB = (f(C) + f(B)) * (B - C) / 2;
		double sACB = sAC + sCB;
		if (abs(sAB - sACB) >= eps * abs(sACB)) {
			SumRecord tmpRecordAC = { A, C, sAC };
			localStack.push(tmpRecordAC);
			A = C;
			sAB = sCB;
		} else {
			J += sACB;
			if (localStack.empty())
				break;
			else {
				tmpRecordAB = localStack.top();
				A = tmpRecordAB.A;
				B = tmpRecordAB.B;
				sAB = tmpRecordAB.s;
				localStack.pop();
			}
		}
	}
	return J;
}

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

Оценить аналитически последовательную сложность представленного алгоритма невозможно, поскольку нельзя указать количество разбиений отрезка интегрирования для каждой отдельно взятой функции, оно зависит от свойств конкретно взятой функции.

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

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

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

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

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

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

Входные данные:

• функция [math]y = f(x)[/math];
• точность [math]\epsilon[/math];
• границы отрезка интегрирования [math]a[/math], [math]b[/math].

Выходные данные:

• число [math]I[/math], равное [math]\int_a^b{f(x)dx}[/math] с точностью [math]\epsilon[/math].

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

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

2 ЧАСТЬ. Программная реализация алгоритма

2.7 Существующие реализации алгоритма

Существует реализация в библиотеке ALGLIB
Есть реализация с использованием OpenMP от университета RWTH AACHEN UNIVERSITY
Кроме того, схема реализации алгоритма представлена Якобовским М. В.