Адаптивная квадратура Гаусса-Кронрода¶
Адаптивная квадратура Гаусса-Кронрода — это численный метод приближённого вычисления определённых интегралов, сочетающий квадратурные формулы Гаусса и Кронрода с автоматическим выбором шага интегрирования (адаптацией) для достижения заданной точности. Метод широко применяется в вычислительной математике, физике и инженерных расчётах благодаря высокой точности и возможности контроля погрешности.
¶История
Квадратура Гаусса была разработана Карлом Фридрихом Гауссом в начале XIX века. Она основана на выборе узлов интегрирования как корней ортогональных многочленов (например, многочленов Лежандра) и обеспечивает максимальную алгебраическую точность для заданного числа узлов. Однако недостатком метода является отсутствие простой оценки погрешности, так как для её вычисления требуется дополнительное увеличение числа узлов, что приводит к пересчёту всех значений подынтегральной функции.
В 1964 году Александр Кронрод предложил модификацию квадратуры Гаусса, добавив к узлам Гаусса дополнительные узлы, расположенные между ними. Это позволило получить квадратурную формулу с более высокой точностью и одновременно оценить погрешность интегрирования. В 1970-х годах идея Кронрода была объединена с адаптивными алгоритмами, что привело к созданию эффективных методов численного интегрирования, реализованных в библиотеках численного анализа (например, QUADPACK).
¶Основные принципы
¶Квадратура Гаусса-Лежандра
Квадратура Гаусса-Лежандра для интеграла $\int_{-1}^{1} f(x) \, dx$ использует $n$ узлов $x_i$ (корни многочлена Лежандра $P_n(x)$) и веса $w_i$, определяемые по формуле: \[ \int_{-1}^{1} f(x) \, dx \approx \sum_{i=1}^{n} w_i f(x_i). \] Эта формула точна для многочленов степени не выше $2n-1$. Для интеграла на произвольном отрезке $[a, b]$ выполняется линейное преобразование.
¶Квадратура Кронрода
Кронрод предложил добавить к $n$ узлам Гаусса ещё $n+1$ узел, расположенный между ними, включая концы отрезка. Общее число узлов становится $2n+1$. Веса для новых узлов вычисляются так, чтобы формула была точна для многочленов степени не выше $3n+1$ (при $n$ нечётном) или $3n$ (при $n$ чётном). Например, для $n=7$ (7 узлов Гаусса) используется 15 узлов Кронрода (7+8). Погрешность оценивается как разность между результатами квадратур Гаусса и Кронрода: \[ E \approx \left| I_{\text{Кронрод}} - I_{\text{Гаусс}} \right|. \]
¶Адаптивный алгоритм
Адаптивная квадратура Гаусса-Кронрода разбивает отрезок интегрирования на подотрезки, на каждом из которых применяется пара квадратур (Гаусса и Кронрода). Если оценка погрешности на текущем подотрезке превышает заданный допуск (абсолютный или относительный), отрезок делится пополам, и процедура повторяется рекурсивно. Алгоритм продолжается до тех пор, пока на всех подотрезках не будет достигнута требуемая точность. Это позволяет эффективно обрабатывать функции с резкими изменениями, осцилляциями или особенностями.
¶Реализация
¶Алгоритм (псевдокод)
`` function adaptive_gauss_kronrod(f, a, b, tol, max_depth): I_gauss, I_kronrod = compute_gauss_kronrod(f, a, b) err = abs(I_kronrod - I_gauss) if err < tol or depth >= max_depth: return I_kronrod else: mid = (a + b) / 2 return adaptive_gauss_kronrod(f, a, mid, tol/2, depth+1) + adaptive_gauss_kronrod(f, mid, b, tol/2, depth+1) ``
¶Популярные реализации
- QUADPACK (Fortran) — библиотека, разработанная в 1970-х годах, включает подпрограмму
dqag(адаптивная квадратура Гаусса-Кронрода). Она используется в GNU Scientific Library (GSL) и SciPy. - SciPy (Python) — функция
scipy.integrate.quadоснована на QUADPACK и использует адаптивную квадратуру Гаусса-Кронрода с 15, 21, 31, 41, 51 или 61 узлом. - MATLAB — функция
integral(с 2012 года) реализует адаптивную квадратуру Гаусса-Кронрода с автоматическим выбором числа узлов. - GSL (C/C++) — функция
gsl_integration_qagпредоставляет адаптивный алгоритм с несколькими правилами Кронрода.
¶Примеры применения
¶Пример 1: Гладкая функция
Вычисление интеграла $\int_{0}^{1} e^{-x^2} \, dx$:
- Точное значение: $\frac{\sqrt{\pi}}{2} \text{erf}(1) \approx 0.746824132812427$.
- Адаптивная квадратура Гаусса-Кронрода (15 узлов) даёт результат с точностью $10^{-12}$ за 1-2 итерации, так как функция гладкая.
¶Пример 2: Функция с особенностью
Вычисление $\int_{0}^{1} \frac{1}{\sqrt{x}} \, dx = 2$:
- Подынтегральная функция имеет особенность в нуле. Адаптивный алгоритм автоматически сгущает узлы вблизи $x=0$, достигая точности $10^{-10}$ за 10-15 делений.
¶Пример 3: Осциллирующая функция
Вычисление $\int_{0}^{10} \sin(100x) \, dx$:
- Высокая частота осцилляций требует большого числа узлов. Адаптивная квадратура Гаусса-Кронрода эффективно разбивает отрезок на мелкие подотрезки, где функция ведёт себя как синусоида.
¶Преимущества и недостатки
¶Преимущества
- Высокая точность: квадратура Кронрода обеспечивает алгебраическую точность до $3n+1$ степени, что превосходит обычную квадратуру Гаусса.
- Автоматический контроль погрешности: оценка погрешности по разности двух квадратур не требует дополнительных вычислений.
- Адаптивность: алгоритм автоматически сгущает узлы в областях с большими градиентами, что экономит вычислительные ресурсы.
- Робастность: метод хорошо работает для большинства гладких и кусочно-гладких функций.
¶Недостатки
- Вычислительная сложность: на каждом шаге требуется вычисление подынтегральной функции во всех узлах (до $2n+1$ точек), что может быть затратно для сложных функций.
- Чувствительность к разрывам: для функций с разрывами первого рода или сингулярностями алгоритм может потребовать очень мелкого разбиения или не сойтись.
- Ограничение на размерность: метод применим только к одномерным интегралам; для многомерных интегралов требуются другие подходы (например, кубатурные формулы или методы Монте-Карло).
¶Сравнение с другими методами
| Метод | Точность | Адаптивность | Оценка погрешности | Сложность |
|---|---|---|---|---|
| Квадратура Гаусса | Высокая (для гладких функций) | Нет | Нет | Низкая |
| Квадратура Кронрода | Очень высокая | Нет | Да | Средняя |
| Адаптивная квадратура Гаусса-Кронрода | Очень высокая | Да | Да | Высокая |
| Метод Симпсона | Средняя | Да | Да | Низкая |
| Метод Монте-Карло | Низкая (сходимость $1/\sqrt{N}$) | Нет | Статистическая | Высокая |
¶Применение
- Физика: вычисление потенциалов, интегралов перекрытия в квантовой механике, интегралов от функций Грина.
- Инженерия: расчёт моментов инерции, центров масс, электрических полей.
- Статистика: численное интегрирование плотностей вероятностей, вычисление моментов распределений.
- Финансовая математика: оценка стоимости опционов (например, метод Блэка-Шоулза требует численного интегрирования).
- Численное решение дифференциальных уравнений: метод Галёркина и метод конечных элементов часто требуют вычисления интегралов от базисных функций.
¶Интересные факты
- Квадратура Кронрода была разработана для улучшения квадратуры Гаусса, но сам Кронрод не использовал адаптивный алгоритм; адаптивность была добавлена позже в библиотеке QUADPACK.
- В библиотеке QUADPACK реализовано несколько правил Кронрода с разным числом узлов: 15, 21, 31, 41, 51 и 61. Обычно используется 15 узлов (7 Гаусса + 8 Кронрода), так как это даёт хороший баланс между точностью и скоростью.
- Адаптивная квадратура Гаусса-Кронрода является стандартным методом для одномерного интегрирования в большинстве научных вычислительных пакетов (SciPy, MATLAB, GSL, NAG).
- Существуют модификации метода для работы с функциями, заданными таблично, или для интегралов с бесконечными пределами (например, преобразование к конечному отрезку).
¶Источники
- Кронрод, А. С. (1964). «О приближённом вычислении определённых интегралов». Доклады АН СССР, 154(2), 283–286.
- Piessens, R., de Doncker-Kapenga, E., Überhuber, C. W., & Kahaner, D. K. (1983). QUADPACK: A Subroutine Package for Automatic Integration. Springer-Verlag.
- Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). Cambridge University Press.
- Документация SciPy:
scipy.integrate.quad(https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html).
BFOmetr — база данных и аналитика по компаниям России.
На главную BFOmetr →
