Квадратурные (кубатурные) методы численного интегрирования по отрезку (многомерному кубу)

Материал из Алговики
Перейти к навигации Перейти к поиску

1 Описание свойств и структуры алгоритма

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

Численное интегрирование — приближенное вычисление значения определенного интеграла. Допустим, нам дана функция [math]f(x)[/math], определенная на отрезке или многомерном кубе, и возможность получать ее численное значение в любой из точек области определения за фиксированное время. Задача — вычислить определенный интеграл данной функции по заданной области и дать оценку погрешности вычисления.

Исторически численное интегрирование также называлось численной квадратурой, а многомерное — кубатурой. Большая советская энциклопедия определяет слово «квадратура» следующим способом: (лат. quadratura — придание квадратной формы), 1) число квадратных единиц в площади данной фигуры. 2) Построение квадрата, равновеликого данной фигуре. 3) Вычисление площади или интеграла.

В данной статье рассматривается семейство методов численного интегрирования, объединенные общим принципом: область интегрирования разбивается по каждой из осей на равное количество частей. В каждой из полученных маленьких областей интеграл приближается простой функцией (например, линейной), значение которой вычисляется явно (путем вычисления подынтегрального выражения в одной или нескольких точках). Ввиду линейности оператора интегрирования по областям полученные значения суммируются и представляют собой результат интегрирования.

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

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

Все перечисленные выше методы, а также подобные им могут быть объединены следующим образом: область интегрирования содержит [math]M[/math] точек, в каждой из которых необходимо вычислить подынтегральное выражение и помножить ее на весовой коэффициент, затем данные числа складываются, и полученная сумма и представляет собой результат интегрирования. Погрешность интегрирования оценивается ровно таким же образом.

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

Исходные данные: функция [math]f[/math] от [math]n[/math] переменных, границы области интегрирования, представляющей собой многомерный куб (по каждой из осей интегрирования заданы минимальное и максимальное значение координаты), максимальное число [math]M[/math] точек, в которых можно вычислить функцию [math]f[/math].

Вычисляемые данные: число, представляющее собой приближенное значение интеграла и, если это позволяет выбранный метод, число, представляющее собой оценку погрешности. По количеству точек [math]M[/math] и методу интегрирования определяются размеры прямоугольных параллелепипедов, на которые разбивается область интегрирования. В каждом из параллелепипедов в зависимости от выбора метода определяются свои точки для вычисления функции [math]f[/math] и весовые коэффициенты.

Наибольшее разнообразие выбора методов достигается в одномерном случае.

Метод прямоугольников: функция для заданного отрезка приближается константной функцией, тем самым, интеграл приближается значением функции в одной точке. Обычно в качестве этой точки используется граница отрезка или, лучше, его середина.

[math]\int\limits_a^b f(x)\,dx \approx f((a+b)/2) * (b-a)[/math]

Метод трапеций: [math]f(x)[/math] приближается линейной функцией, тем самым, интеграл приближается усредненным значением по краям отрезка:

[math]\int\limits_a^b f(x)\,dx \approx (f(a) + f(b)) * (b-a)/2[/math]

Метод парабол (метод Симпсона): [math]f(x)[/math] приближается параболой, пересекающей график функции на краях отрезка и в его середине:

[math]\int\limits_a^b f(x)\,dx \approx \frac{b-a}{6}\left(f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right)[/math]

Метод Гаусса: [math]f(x)[/math] приближается прямой как в методе трапеций, однако в качестве точек выбираются не концы отрезка, что позволяет увеличить точность:

[math]\int\limits_a^b f(x)\,dx \approx \frac{b-a}{2}\left(f\left(\frac{a+b}{2} - \frac{b-a}{2\sqrt{3}} \right)+f\left(\frac{a+b}{2} + \frac{b-a}{2\sqrt{3}} \right) \right)[/math]

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

Так или иначе, на подготовительной стадии определяется M наборов координат [math]{x_{1,i}, … x_{n,i}}[/math], где [math]i[/math] меняется от [math]1[/math] до [math]M[/math], в которых будет вычисляться функция [math]f[/math] и набор весовых коэффициентов [math]с_i[/math], на которые будет умножаться значение функции, после чего управление передается вычислительному ядру алгоритма.

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

Вычислительное ядро последовательной версии любого из квадратурных методов численного интегрирования представляет собой вычисление суммы

[math]\sum \limits_1^M c_i f(x_{1,i}, … x_{n,i}) [/math]

Поскольку время, требуемое на операции сложения и умножения в реальных примерах представляют собой о малое от времени, затрачиваемого на вычисление функции [math]f[/math], для производительности становится совершенно неважно, каким способом осуществлять сложение и умножение.

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

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

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

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

1.6 Последовательная сложность алгоритма

Время, затрачиваемое на вычисление функции растет примерно как [math]t*M[/math], где [math]t[/math] — время, требуемое на вычислении функции [math]f[/math], а [math]M[/math] — количество точек, в которых вычисляется функция [math]f[/math]. Тем самым, алгоритм является линейным относительно количества точек интегрирования.

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

Опишем граф алгоритма как аналитически, так и в виде рисунка.

Граф алгоритма состоит из трех групп вершин. Вершины на верхнем ярусе берут в качестве входных данных координаты и вычисляют по ним значение функции [math]f[/math]. Вершины на втором ярусе умножают полученное значение на соответствующий весовой коэффициент. Вершины на третьем ярусе поочередно складывают полученные значения; результат относится к выходным данным.

НА ЧЕМ ДЕЛАЕТСЯЯ РИСУНОК???????

1.8 Описание ресурса параллелизма алгоритма

Для вычисления интеграла функции методом кубатур требуется последовательно выполнить следующие ярусы: - Вычисление исходной функции [math]f[/math], производящееся независимо для различных точек и выполняющееся за близкое к константному время [math]t[/math] — основная часть алгоритма. - Умножение полученных результатов на соответствующие весовые коэффициенты - Сложение полученных результатов.

В общем случае время, затрачиваемое на сложение или умножение пренебрежимо мало по сравнением со временем необходимым на вычисление функции [math]f[/math], а все вычисления функции [math]f[/math] могут производиться на различных узлах, тем самым, задача вычисления интеграла методом кубатур идеально поддается параллельной реализации. Однако при значительном увеличении количества узлов (в случае, когда оно становится сравнимым с количеством точек, в которых вычисляется функция), ресурс параллелизма может быть ограничен накладными расходами,такими как передача исходного выражения на каждый из вычислительных узлов.

Таким образом, до тех пор, пока количество вычислительных узлов не превышает количества точек, в которых должно быть произведено интегрирование, алгоритм выполняется за около константное время. Однако, в реальных вычислениях обычно требуется значительно более высокое количество точек интегрирования (измеряемое миллионами), чем количество узлов, которые можно задействовать, и в таком случае сложность алгоритма является линейной.