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

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

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

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