Main menu

Численное решение интегральных уравнений: методы Фредгольма и Вольтерры

Когда неизвестная функция скрыта под знаком интеграла

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

Интегральные уравнения принято классифицировать по двум признакам. Если пределы интегрирования фиксированы (постоянны), то это уравнение Фредгольма. Если верхний предел является переменной интегрирования (то есть интеграл зависит от текущего положения), то это уравнение Вольтерры, описывающее процессы эволюции с памятью. Если неизвестная функция находится только под знаком интеграла, это уравнение первого рода; если она присутствует еще и вне интеграла — это уравнение второго рода. Как и в случае с ОДУ, аналитически решить интегральные уравнения удается лишь в редких частных случаях, что делает применение численных методов абсолютно необходимым.

Метод квадратур: сведение к СЛАУ

Самым прямым и универсальным подходом к решению интегральных уравнений Фредгольма второго рода является метод квадратур. Идея этого метода базируется на классическом численном интегрировании, которое мы обсуждали ранее. Интеграл, содержащий неизвестную функцию, приближенно заменяется конечной суммой — квадратурной формулой (например, по правилу трапеций или Симпсона) на дискретной сетке узлов.

В результате такой дискретизации исходное непрерывное уравнение для бесконечного числа точек превращается в систему линейных алгебраических уравнений (СЛАУ) относительно неизвестных значений функции в узлах сетки. Матрица этой СЛАУ, как правило, является плотной (полностью заполненной ненулевыми элементами), так как интегральный оператор связывает все точки области (глобальное взаимодействие). Решив эту плотную СЛАУ стандартным методом Гаусса или QR-разложением, мы получаем приближенные значения искомой функции в узлах. Точность метода напрямую зависит от плотности сетки и порядка используемой квадратурной формулы.

Метод последовательных приближений (итерации Пикара)

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

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

Подробнее

Адаптивные методы интегрирования ОДУ: алгоритмы Рунге-Кутты-Фельберга

Компромисс между скоростью и точностью: проблема постоянного шага

Решение обыкновенных дифференциальных уравнений (ОДУ) методами Рунге-Кутты, которые мы рассматривали ранее, часто проводится с постоянным шагом интегрирования. Однако в реальных физических системах характер процессов может резко меняться с течением времени. Представьте космический аппарат, который летит к Луне. Большую часть пути гравитационные поля меняются очень плавно, и траекторию можно рассчитывать с огромным шагом (например, раз в час). Но при подлете к орбите Луны и выполнении маневров гравитационные силы начинают меняться стремительно. Если оставить прежний большой шаг, расчетная траектория грубо отклонится от реальной. Если же с самого начала установить микроскопический шаг «с запасом» для прохождения сложных участков, то расчет простого участка пути займет неоправданно много времени процессора.

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

Вложенные методы и оценка локальной погрешности

Блестящее решение этой проблемы предложил немецкий математик Эрвин Фельберг в 1969 году. Он разработал концепцию вложенных методов Рунге-Кутты. Идея заключается в том, чтобы на каждом шаге вычислять новое приближение функции сразу двумя методами разного порядка точности (например, четвертого и пятого порядка).

Если бы мы делали это «в лоб», объем вычислений возрос бы в два раза. Гениальность идеи Фельберга (и последующих авторов, таких как Дорман и Принс) состоит в том, что коэффициенты для методов 4-го и 5-го порядков подбираются специальным образом. Они используют одни и те же промежуточные вычисления (те самые коэффициенты наклона касательных внутри шага). То есть, вычислив набор наклонов для метода 5-го порядка, мы можем путем простого перемножения с другими весами получить и результат для 4-го порядка практически бесплатно! Разность между результатами 5-го и 4-го порядка служит превосходной оценкой локальной ошибки (погрешности усечения) на данном шаге.

Алгоритм управления шагом (Step Size Control)

Имея оценку локальной ошибки, программа-решатель включает встроенный ПИД-регулятор шага. Пользователь заранее задает желаемую допустимую погрешность (tolerance). Алгоритм сравнивает вычисленную оценку ошибки с заданным допуском.

Если ошибка оказалась больше допуска, это означает, что шаг был слишком велик. Алгоритм отбрасывает вычисленное значение, уменьшает шаг по специальной формуле (обычно умножая на коэффициент около 0.8 от расчетного оптимального, чтобы иметь запас надежности) и пересчитывает этот участок заново. Если же ошибка оказалась значительно меньше допуска, алгоритм принимает шаг как успешный, но для следующего шага автоматически увеличивает размер интервала, ускоряя симуляцию. Метод Рунге-Кутты-Фельберга (RKF45) и его более современная модификация, метод Дормана-Принса (DOPRI5), стали мировым стандартом адаптивного интегрирования и используются по умолчанию во встроенных функциях MATLAB (ode45) и библиотеках SciPy в Python.

Подробнее

Методы многомерной оптимизации: градиентный спуск и метод Ньютона

Поиск наилучшего решения в мире множества вариантов

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

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

Градиентный спуск: движение по наискорейшему склону

Самым известным и концептуально простым алгоритмом является метод градиентного спуска (Gradient Descent), предложенный еще Огюстеном Луи Коши. Градиент — это вектор, указывающий направление наискорейшего возрастания функции. Следовательно, антиградиент (вектор, взятый со знаком минус) указывает направление наискорейшего убывания.

Алгоритм стартует из случайной начальной точки. На каждом шаге мы вычисляем вектор градиента в текущей точке и делаем небольшой шаг в противоположном направлении. Размер этого шага определяется параметром, который называется «скоростью обучения» (learning rate). Главная проблема классического градиентного спуска — сильная зависимость от этого параметра. Если шаг слишком велик, алгоритм начнет перепрыгивать через минимум и расходиться; если слишком мал — сходимость растянется на долгие часы. Кроме того, в так называемых «овражных» функциях градиентный спуск начинает совершать зигзагообразные колебания поперек оврага, крайне медленно продвигаясь к его дну.

Метод Ньютона и учет кривизны пространства

Чтобы ускорить сходимость и избежать зигзагов в оврагах, необходимо использовать информацию не только о наклоне функции (первой производной), но и о ее кривизне. Эту информацию предоставляет матрица вторых производных, называемая матрицей Гессе (или гессианом). Использование гессиана лежит в основе многомерного метода Ньютона.

Метод Ньютона строит на каждом шаге квадратичную аппроксимацию целевой функции (в виде параболоида) и мгновенно прыгает в вершину этого параболоида. В окрестности минимума метод Ньютона обладает фантастической квадратичной скоростью сходимости. Однако вычисление, хранение и, что самое главное, обращение матрицы Гессе на каждом шаге (требующее решения огромной СЛАУ) делает классический метод Ньютона абсолютно неприменимым для задач с большим количеством параметров (например, при обучении нейросетей). Поэтому на практике используют квазиньютоновские методы (такие как популярный BFGS) или метод сопряженных градиентов, которые хитроумными способами оценивают кривизну пространства, не вычисляя сам гессиан в явном виде.

Подробнее

Соц. сети