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

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

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_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.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.8 Ресурс параллелизма алгоритма

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

пока в глобальном стеке есть отрезки:

1. взять один отрезок из глобального стека;
2. выполнить алгоритм локального стека, выполняя при записи в локальный стек следующие операции:
   1. если

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 Существующие реализации алгоритма