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

Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки
Последовательный алгоритм
Объём входных данных	[math]O(1)[/math]
Объём выходных данных	[math]O(1)[/math]

Основные авторы описания: Соколов К. В. группа 604

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

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

Основаная идея метода заключается в разбиении области интегрирования на небольшие отрезки, в каждом из которых подынтегральная функция заменяется на константу (метод прямоугольников), линейную фунцкию(метод трапеций), полином второй степени(метод Симпсона). Адаптивность сетки означает, что если на одном из сегментов разбиения приближенное значение интеграла, вычисленное на сетке мелкости [math]h[/math] значительно(более чем на [math]\varepsilon[/math]) отличается от приближенного значения интеграла, вычисленного на этом же сегменте, но на сетке мелкости [math]\frac{h}{2}[/math], в качестве приближенного значения интеграла мы должны выбрать то, что вычислено на более мелкой сетке и рекурсивно проверить, вычисленно ли оно с достаточной точностью.

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

На вход мы получаем функцию [math]f(x)[/math], левую границу отрезка [math]A[/math], правую границу отрезка интегрирования [math]B[/math], точность [math]\varepsilon[/math]. Требуется вычислить значение определенного интеграла функции [math]f(x)[/math]:

[math]\int_A^B f(x) dx[/math]

Разобьем область интегрирования [math][A; B][/math] на некоторое число равных отрезков [math]N[/math]. Получим последовательность узлов разбиения [math][A = x_0 \lt x_1 \lt .. \lt x_{N} = B][/math]. Тогда можно вычислить приближенное значение интеграла:

[math]I(f(x)) = \sum_{i=0}^{N-1} J(x_i, x_{i + 1})[/math], где [math]J(x_i, x_{i + 1}) = S(x_i, x_{i + 1})[/math] если [math]S(x_i, x_{i + 1}) - S(x_i, x_{i + \frac{1}{2}}) - S(x_{i + \frac{1}{2}}, x_{i + 1}) \lt \varepsilon[/math],

иначе [math]J(x_i, x_{i + 1}) = J(x_i, x_{i + \frac{1}{2}}) + J(x_{i + \frac{1}{2}}, x_{i + 1})[/math]

[math]S(x_i, x_j) = \frac{f(x_i) + f(x_j)}{2} * (x_j - x_i)[/math]

[math]x_{i + \frac{1}{2}} = \frac{x_i + x_{i + 1}}{2}[/math]

Вычисленное предложенным образом [math] I(f(x))[/math] будет приближенным значением интеграла, вычисленным с точностью [math]\varepsilon[/math]

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

Вычислительным ядром алгоритма является вычисление приближенного значения интеграла на [math]i[/math]-ом сегменте разбиения [math]J(i, i + 1)[/math]. Важным моментом является сложность вычисления подынтегральной функции [math]f(x)[/math], если это сложность достаточно велика, алгоритм можно реализовать с минимальным количеством пересчетов значений [math]f(x)[/math]

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

Для реализации алгоритма предложен метод локального стека без ограничения глубины рекурсии.

Суть метода заключается в том, что все отрезки разбиения хранятся в стеке и в случае если на текущем отрезке достигнута необходимая точность интегрирования, то приближенное значение на текущем отрезке добавляется к общей сумме и берется следующий отрезок из стека(если стек пуст, то ответ получен, выполнение прекращается). Иначе отрезок разбивается на две части, правую часть кладем в стек, на левой вычисляем приближенное значение интеграла.

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

Реализация описана в книге Введем стек с тремя операциями:

   sp=0  //    указатель вершины стека
   struct {
         A, B, fA, fB, s
   } stk[1000] 

1. PUT_INTO_STACK - положить элемент в стек

   #define PUT_INTO_STACK(A, B, fA, fB, s)
   {
       stk[sp].A=A
       stk[sp].B=B
       stk[sp].fA=fA
       stk[sp].fB=fB
       stk[sp].s=s
       sp++
   } 

2. GET_FROM_STACK - взять верхний элемент из стека

   #define GET_FROM_STACK(A, B, fA, fB, s)
   {
       sp--
       A=stk[sp].A
       B=stk[sp].B
       fA=stk[sp].fA
       fB=stk[sp].fB
       s=stk[sp].s
   }

3. STACK_IS_NOT_FREE - проверить стек на наличие в нем элементов

   #define STACK_IS_NOT_FREE (sp > 0)

Непосредственное описание самого алгоритма:

   LocalStack (A, B)
   {
       // В переменной J храним ответ
       J=0
       // Вычисляем значение функции на границах интегрируемого отрезка
       fA = f(A)
       fB = f(B)
       // Вычисляем площадь с помощью метода трапеций
       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
           sACB = sAC + sCB
           if (|sAB - sACB| ≥ ε|sACB|) {  // Если заданная точность не достигнута
               PUT_INTO_STACK(A, C, fA, fC, sAC)  // то один из отрезков кладем в стек
               A = C  // и решаем задачу для оставшегося отрезка
               fA = fC
               sAB = sCB
           } else  // Если заданная точность достигнута
           {
               J += sACB  // то к общей сумме добавляем полученное число
               if (STACK_IS_NOT_FREE)  // если стек пуст, то цикл завершается - ответ получен
                   break 
               GET_FROM_STACK(A, B, fA, fB, sAB)  // иначе берем из стека новый отрезок
           } 
       }
       return J 
   }

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

Глубина рекурсии при интегрировании существенно зависит от вида входной функции, однако в случае ограниченной глубины рекурсии(что приведет к потере гарантированной точности) сложность алгоритма можно оценить как [math]O(N)[/math]

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

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

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

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

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

Вход: Границы сегмента интегрирования, точность (3 числа типа double) [math]\varepsilon[/math]

Выход: приближенное значение определенного интеграла (Число типа double)

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

1. Алгоритм балансирует нагрузку между процессами, в отличии от, например, метода геометрического параллелизма

2. Алгоритм можно изменить в зависимости от того, что дороже - пересылать значения функции между процессами или пересчитывать подынтегральную функцию снова

3. Алгоритм устойчив

4. Алгоритм детерминированный

5. В идеальном случае(линейная функция) алгоритм даст точное значение интеграла за один шаг, отношение сложности параллельного и последовательного алгоритмов составит n.

6. Вычислительная сложность алгоритма существенно зависит от вида подынтегральной функции, для отрезка длины [math]L[/math] ее можно оценить как [math]\frac{L}{N}k[/math], где [math]k[/math] - количество операций для вычисления "самого сложного" (требующего наибольшего количества делений) отрезка из всей области, [math]N[/math] - число процессоров, выполняющих интегрирование

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

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

Реализация с использованием OpenMP [1]

3 Литература

<references \>