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

Основные авторы описания: Можарова В. А.

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

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

Основная идея большинства методов численного интегрирования состоит в замене подынтегральной функции на более простую, интеграл от которой легко вычисляется аналитически. При замене подынтегральной функции на полином нулевой, первой и второй степени получаются соответственно методы прямоугольников, трапеций и парабол (Симпсона).

Чтобы посчитать определенный интеграл, участок интегрирования делится на отрезки, на которых подсчитывается значение новой (упрощенной) функции, после чего значения суммируются на всех отрезках. Но чтобы посчитать интеграл с определенной точностью, мы не можем заранее узнать, какие размеры отрезков, нужно использовать, т. к. это зависит от самой функции. Для этого используют адаптивное разбиение участка интегрирования: размер конкретного отрезка определяется во время подсчета значения функции на нем самом. Если заданная точность не была достигнута на данном отрезке, он рекурсивно разбивается на несколько новых и функция подсчитывается уже на них. Таким образом, выигрыш от применения рассмотренного алгоритма рекурсивного разбиения достигается за счет возможности использования сеток с разным числом узлов на разных участках интервала интегрирования.

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

Пусть у нас есть интеграл, который надо вычислить с точностью ε:

Предположим, что в результате адаптации сетки для данной функции, на отрезке [A, B] получилась сетка, содержащая n + 1 узел. Тогда согласно методу трапеций, можно численно найти определенный интеграл от функции f(x) на отрезке [A, B]:

Будем полагать значение J найденным с точностью ε, если выполнено условие:

где n₁ и n₂ - количество узлов, образующих две разные сетки. Пользуясь формулой

рекурсивно делим отрезки интегрирования до достижения на них заданной точности.

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

Основное время работы алгоритма приходится на вычисление площадей на отрезках, получившихся в результате рекурсивного разбиения:

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

Для хранения отрезков интегрирования в данном алгоритме понадобится стек .

   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)

1.5 Схема реализации алгоритма

   LocalStack (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.7 Информационный граф

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

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

На вход алгоритму подается подынтегральная функция, пределы интегрирования и значение точности, с которой нужно посчитать заданный интеграл. На выходе: число – значение определенного интеграла.

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

1. Алгоритм адаптируется под сложность функции

2. Для каждого узла сетки значение функции высчитывается один раз

3. В параллельной версии алгоритма не требуется управляющий процесс

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

5. Степень исхода вершины информационного графа равна 1, т. к. результат работы каждого процесса в дальнейшем используется только 1 раз.

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

7. Данный алгоритм является устойчивым

3 Литература

<references \>