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

Вычисление определенного интеграла с использованием адаптивно сгущающейся сетки
Последовательный алгоритм
Последовательная сложность	[math] O\left( n + 2^{log 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 Общее описание алгоритма

---

Пусть требуется вычислить определённый интеграл [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.3 Вычислительное ядро алгоритма

---

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

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

---

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

   int sp = 0;       //указатель вершины стека

   struct
   {
       double a;
       double b;
       double fa;
       double fb;
       double sab;
   } stk[1000];      //массив структур, в которых хранятся отложенные значения. 1000 может быть заменена на необходимое число

   //макроопределения доступа к данным стека
   #define STACK_IS_NOT_FREE ( sp > 0 )

   #define PUT_INTO_STACK( a, b, fa, fb, sab )
   {
       stk[sp].a = a;
       stk[sp].b = b;
       stk[sp].fa = fa;
       stk[sp].fb = fb;
       stk[sp].sab = sab;
       ++sp;
   }

   #define GET_FROM_STACK( a, b, fa, fb, sab )
   {
       --sp;
       a = stk[sp].a;
       b = stk[sp].b;
       fa = stk[sp].fa;
       fb = stk[sp].fb;
       sab = stk[sp].sab;
   }

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

---

Будем считать, что описана функция double f( double x ), возвращающая значение f(x) интегрируемой функции f в точке x.

   double f( double x );

   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 c = 0;        //середина отрезка интегрирования
       double sab = 0;      //
       double sac = 0;      //
       double scb = 0;      //
       double sabc = 0;     //переменные для подсчёта текущих сумм
       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 - sacb ) >= [math]\epsilon[/math] * abs ( sacb ) ) //[math]\epsilon[/math] - заданная точность
           {
               PUT_INTO_STACK( a, c, fa, fc, sac );
               a = c; 
               fa = fc;
               sab = scb;
           }
           else
           {
               I += sacb;
               if( STACK_IS_NOT_FREE )
               {
                   break;
               }
               GET_FROM_STACK( a, b, fa, fb, sab );
           }
       }
       return 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) [/math].

При этом из всего вышесказанного сразу видно, что при достаточном количестве вычислительных ресурсов распараллеливания высота ярусно-параллельной формы параллельной программы будет равна глубине [math] k \leq log(n), [/math], а ширина будет не больше, чем количество отрезков [math] n. [/math]

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].

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

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