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

Авторы: Устинов Г.А. Лиджиев Д.В.

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 Макроструктура алгоритма

Метод локального стека

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

   // макроопределения доступа к данным стека

   #define STACK_IS_NOT_FREE (sp>0)

   #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++
   } 

   #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
   }

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

   IntTrap (A, B)
   {
       J=0
       fA = f(A)
       fB = f(B)
       sAB = (fA + fB) * (B - A) / 2
       while (1)
       {
           C = (A + B) / 2
           fc = fun(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 Последовательная сложность алгоритма

.

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

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

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

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. В параллельной версии алгоритма на основе не требуется управляющий процесс. Распределение отрезков и финальное сложение можно провести на нулевом процессе.

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

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

Существует реализация в библиотеке ALGLIB.

3 ЧАСТЬ. Литература

Якобовский М.В., Кулькова Е.Ю. Решение задач на многопроцессорных вычислительных системах с разделяемой памятью. - М.: СТАНКИН, 2004. – 30 c.