Участник:Бобцов Борис/Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки
| Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки | |
| Последовательный алгоритм | |
| Последовательная сложность | [math] O\left( n \right) [/math], где [math] n [/math] - количество подотрезков, на которые разбился исходный отрезок |
| Объём входных данных | [math]1[/math] функция [math]1[/math] переменной, [math]3[/math] числа (границы отрезка интегрирования, точность) |
| Объём выходных данных | [math]1[/math] число |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math] log(n) [/math] |
| Ширина ярусно-параллельной формы | [math] n [/math] |
Автор статьи: Бобцов Борис (группа 619)
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Пусть требуется вычислить определённый интеграл [math]I = \int_a^b{f(x)dx}[/math] от непрерывной функции [math]f(x)[/math]. Если может быть найдена первообразная [math]F(x)[/math] подынтегральной функции, то по формуле Ньютона-Лейбница:
[math]\int_a^b{f(x)dx} = F(b) - F(a). [/math]
Если же первообразная не может быть найдена, то для вычисления интеграла прибегают к приближённым формулам, точность которых может быть сделана сколь угодно большой.
Приближённые методы вычисления определённого интеграла в большинстве случаев основаны на том, что определённый интеграл численно равен площади криволинейной трапеции, ограниченной кривой [math] y =f(x)[/math], сегментом [math] \left[ a,b \right] [/math] оси [math]Ox[/math] и вертикальными прямыми [math]x = a[/math] и [math]x = b[/math].
Идея данных методов заключается в приближении исходной кривой новой, достаточно близкой к старой, кривой. В качестве этой новой кривой подбирают такую, для которой площадь криволинейной трапеции подсчитывается просто. В зависимости от выбора новой кривой мы получим ту или иную приближённую формулу интегрирования.
Разбивая сегмент интегрирования на несколько достаточно малых частей, на каждом из полученных маленьких сегментов мы можем приблизить интегрируемую функцию прямой. Тем самым мы заменим криволинейную трапецию прямолинейной, формула площади которой проста для вычисления.
1.2 Математическое описание алгоритма
1.2.1 Метод трапеций
Будем рассматривать проблему вычисления значения определенного интеграла с некоторой заданной точностью [math] \epsilon [/math]:
[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]
Будем полагать значение [math] I [/math] найденным с точностью [math] \epsilon [/math], если выполнено условие:
[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.2.2 Метод рекурсивного деления
Метод трапеций обладает следующими недостатками:
- в некоторых точках значение функции будет вычисляться несколько раз;
- на всём интервале интегрирования используется равномерная сетка, тогда как число точек, необходимых для достижения заданной точности, зависит от вида функции и может сильно различаться на разных участках сегмента интегрирования.
Метод рекурсивного деления был придуман, чтобы избежать этих недостатков. Суть метода заключается в рекурсивном вызове себя до момента, пока мы не достигнем нужной точности: [math] I(a,b) = IntTrap(a, b, f(a), f(b)) = IntTrap(a, (b - a) / 2, f(a), f((b - a) / 2)) + IntTrap((b - a) / 2, b, f((b - a) / 2), f(b)) = ... [/math]
Однако можно заметить, что координаты концов обрабатываемых отрезков хранятся в программном стеке процесса [math]IntTrap() [/math], и недоступны программисту. Это является основной сложностью при построении параллельного алгоритма на основе метода рекурсивного деления.
1.2.3 Метод локального стека
Данный алгоритм использует массив данных, доступ к которым осуществляется по принципу стека - первым удаляется последний из добавленных элементов. Теперь у разных процессов есть возможность обрабатывать отрезки независимо друг от друга. В рассматриваемом алгоритме вместо стека можно было использовать и другие способы хранения заданий, например очередь. Стек выбран исключительно из соображений простоты реализации. Алгоритм локального стека можно считать наилучшим из перечисленных последовательных методов для разработки его параллельной версии.
1.3 Вычислительное ядро алгоритма
В схеме реализации будет явно видно, что алгоритм реализуется при помощи цикла, внутри которого происходит подсчет текущих частичных сумм и проверка достигнутой точности. Сложность данного алгоритма зависит от количества проходов по циклу линейно, поэтому можно считать, что тело этого цикла и является вычислительным ядром данного алгоритма.
1.4 Макроструктура алгоритма
Как уже было замечено, алгоритм представляет собой три основные части, выполняющиеся в цикле. Сначала вычисляются частичные суммы на имеющемся отрезке, затем проверяется точность вычислений и, в случае, если точность недостаточна, вычисляются новые отрезки.
1.5 Схема реализации последовательного алгоритма
Будем считать, что описана функция double f( double x ), возвращающая значение f(x) интегрируемой функции f в точке x.
double f( double x );
int sp = 0; //указатель вершины стека
struct
{
double a;
double b;
double fa;
double fb;
double sab;
} stk[1000]; //массив структур, в которых хранятся отложенные значения. 1000 может быть заменена на необходимое число
//операции доступа к данным стека #define STACK_IS_NOT_FREE ( sp > 0 )
void PUT_INTO_STACK( double a, double b, double fa, double fb, double sab )
{
stk[sp].a = a;
stk[sp].b = b;
stk[sp].fa = fa;
stk[sp].fb = fb;
stk[sp].sab = sab;
++sp;
}
void GET_FROM_STACK( double& a, double& b, double& fa, double& fb, double& sab )
{
--sp;
a = stk[sp].a;
b = stk[sp].b;
fa = stk[sp].fa;
fb = stk[sp].fb;
sab = stk[sp].sab;
}
//основная функция
int IntTrap(a, b)
{
double I = 0; //значение интеграла
double fa = f( a ); //значение функции [math]y = f(x)[/math] в точке [math]a[/math]
double fb = f( b ); //значение функции [math]y = f(x)[/math] в точке [math]b[/math]
double fc = 0;
double c = 0; //середина отрезка интегрирования
double sab = 0; //
double sac = 0; //
double scb = 0; //
double sabc = 0; //переменные для подсчёта текущих сумм
sab = ( fa + fb ) * ( b - a ) / 2;
while( 1 )
{
c = ( a + b ) / 2;
fc = f( c );
sac = ( fa + fc ) * ( c - a ) / 2;
scb = ( fc + fb ) * ( b - c ) / 2;
sabc = sac + scb;
if( abs( sab - sabc ) > [math]\epsilon[/math] * abs ( sabc ) ) //[math]\epsilon[/math] - заданная точность
{
PUT_INTO_STACK( a, c, fa, fc, sac );
a = c;
fa = fc;
sab = scb;
}
else
{
I += sabc;
if( STACK_IS_NOT_FREE )
{
break;
}
GET_FROM_STACK( a, b, fa, fb, sab );
}
}
return I;
}
Работа со стеком происходит через операции доступа к данным. В основной функции в цикле происходит вычисление частичных сумм sac и scb, затем проверяется точность if( abs( sab - sacb ) >[math]\epsilon[/math] * abs ( sacb ) ). В случае достижения необходимой точности частичные суммы добавляются к значению интеграла I, если отрезков не осталось, то программа завершает работу, иначе из стека выгружается следующий отрезок. Если точность не достигнута, то правый полуотрезок загружается в стек, а левый полуотрезок устанавливается как текущий для следующей итерации цикла.
Когда в стеке не останется отрезков, это будет означать, что программа рассчитала все отрезки до точки a включительно, и возвращает значение I как приближённое значение интеграла.
1.6 Последовательная сложность алгоритма
Пусть для работы алгоритма было необходимо разбить отрезок [math] \left[a, b\right] [/math] на сегменты [math] \left[a=x_0 , x_1\right], \left[a=x_1 , x_2\right],... \left[a=x_{n-1} , x_n = b\right] [/math]. Так как каждый раз, когда нам приходится вычислять площадь трапеции на более коротком сегменте, мы делим отрезок пополам, получаем, что длина каждого отрезка равна [math] \left( x_i - x_{i - 1} \right) = \frac{b - a}{2^k} [/math] для некоторого натурального [math] k [/math].
Назовём глубиной метода максимальное [math] k [/math] из всех [math] k_i : \left( x_i - x_{i - 1} \right) = \frac{b - a}{2^{k_i}} [/math].
Очевидно, что для каждой функции [math] y [/math] и для каждой точности [math] \epsilon [/math] глубина алгоритма будет разной. То есть можно считать, что [math] k = k(f(x), \epsilon). [/math] Также понятно, что, если [math] k [/math]- глубина алгоритма при фиксированных [math] y [/math] и [math] \epsilon [/math], то максимально возможное количество полученных отрезков, для которых считалась частичная сумма, и количество точек на отрезке [math] \left[a, b\right] [/math] равны соответственно [math] 2^k [/math] и [math] 2^k + 1 [/math] при [math] k_i = k, i = 1,2,...,n [/math]. В этом худшем с точки зрения количества операций случае алгоритм прокрутит цикл while( 1 ) [math] 1 + 2 + 4 + ... + 2^k = (2^{k + 1}-1) [/math] раз.
Если не учитывать при расчёте сложности операции работы с памятью, а только сложения/вычитания и умножения*деления, получим, что за каждый выполненный цикл while( 1 ) алгоритм выполняет 7 сложений и 6 умножений. Итого получаем в худшем случае [math] 7*2^k [/math] сложений и [math] 6*2^k [/math] умножений. Кроме этого, алгоритм посчитает значение [math] f(x) [/math] в [math] n+1 [/math] точках. При этом понятно, что [math] k \leq log(n). [/math]
Получаем, что при разбиении исходного отрезка на [math] n [/math] подотрезков получаем глубину алгоритма [math] k \leq log(n) [/math] и сложность алгоритма [math] O\left( n + 2^{log n} \right) = O\left( n + n \right) = O\left( n \right) [/math].
При этом из всего вышесказанного сразу видно, что при достаточном количестве вычислительных ресурсов распараллеливания высота ярусно-параллельной формы параллельной программы будет равна глубине [math] k \leq log(n), [/math], а ширина будет не больше, чем количество отрезков [math] n. [/math]
1.7 Информационный граф
На рисунке изображён информационный граф алгоритма. Синими прямоугольниками изображены входные данные итерации цикла, внутри подписаны имена данных. Кругами изображены операции вычисления переменных, имена которых указаны внутри кругов. Зелёными овалами изображены прочие операции, выполняемые алгоритмом. Push и Pop - операции добавления и взятия из стека. После операции проверки точности алгоритм переходит на операцию Push, если требуемая точность не достигнута, и на операцию расчёта переменной I - иначе. Проверка стека подразумевает проверку на наличие непросчитанных отрезков. При отсутствии таковых итерация завершается, иначе из стека загружается новый отрезок.
1.8 Ресурс параллелизма алгоритма
Представленный алгоритм обладает хорошим ресурсом параллелизма. В самом деле, пусть в нашем распоряжении есть доступный всем параллельным процессам список отрезков интегрирования, организованный в виде стека. Назовём его глобальным стеком. Предположим, что перед запуском параллельных процессов в глобальный стек помещается единственная запись - координаты концов отрезка интегрирования - в дальнейшем - "отрезок". Тогда возможно реализовать следующую схему алгоритма, выполняемую каждым из параллельных процессов:
Пока в глобальном стеке есть отрезки:
- взять один отрезок из глобального стека;
- выполнить алгоритм локального стека, выполняя при записи в локальный стек следующие операции:
- если в локальном стеке есть несколько отрезков, а в глобальном стеке отрезки отсутствуют, то
- переместить часть отрезков из локального стека в глобальный;
- по исчерпании локального стека добавить полученную частичную сумму к общему значению интеграла.
- если в локальном стеке есть несколько отрезков, а в глобальном стеке отрезки отсутствуют, то
1.9 Входные и выходные данные алгоритма
Входные данные:
- функция [math] y = f(x)[/math];
- точность [math] \epsilon [/math];
- границы отрезка интегрирования [math] a, b.[/math]
Выходные данные:
- число [math]I [/math], приближённо равное [math] \int_a^b{f(x)dx}[/math] с точностью [math] \epsilon [/math].
1.10 Свойства алгоритма
Алгоритм обладает большой вычислительной мощностью. Кроме того, у него имеется большой потенциал для распараллеливания.
Алгоритм не обладает свойством детерминированности, так как в зависимости от специфических свойств графика входных функций может вести себя для них совершенно по-разному.
Соотношение последовательной и параллельной сложности алгоритма высоко и равно приблизительно [math] \frac{n}{log(n)} [/math].
2 Программная реализация алгоритма
2.7 Существующие реализации алгоритма
На данный момент реализация алгоритма локального стека была указана только в книге Якобовского М.В., Кульковой Е.Ю. "Решение задач на многопроцессорных вычислительных системах с разделяемой памятью."
3 Литература
- Шнейдер В. Е. и др. Краткий курс высшей математики. Учеб. пособие для втузов. М., «Высш. школа», 1972. 640 с.
- Якобовский М.В., Кулькова Е.Ю. Решение задач на многопроцессорных вычислительных системах с разделяемой памятью. - М.: СТАНКИН, 2004. – 30 c.
