Main menu

Решение жестких краевых задач: метод ортогональной прогонки Годунова

Вычислительная катастрофа при интегрировании

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

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

Спасение через ортогонализацию

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

Алгоритм интегрирует систему независимых векторов с левого конца отрезка. Как только программа замечает, что векторы начинают «слипаться» из-за влияния жесткости (угол между ними становится слишком острым, а матрица их решений теряет обусловленность), процесс интегрирования останавливается. К текущим векторам применяется классический процесс ортогонализации Грама-Шмидта. Векторы принудительно разворачиваются под прямыми углами друг к другу, а их длины нормируются, предотвращая арифметическое переполнение. При этом алгоритм запоминает матрицы преобразования (матрицы ортогонализации). После «очистки» векторов интегрирование продолжается до следующей опасной точки. Дойдя до правого конца и удовлетворив второе краевое условие, алгоритм совершает обратный ход, используя сохраненные матрицы для восстановления истинного физического решения. Метод Годунова стал золотым стандартом надежности для расчета напряженно-деформированного состояния сложнейших оболочечных конструкций.

Подробнее

Сингулярное разложение матриц (SVD): математика рекомендательных систем и сжатия данных

Универсальный инструмент матричного анализа

В вычислительной линейной алгебре существует множество способов разложения матриц (LU, QR, спектральное разложение), однако Сингулярное Разложение (Singular Value Decomposition, SVD) занимает среди них абсолютно особое место. Если спектральное разложение применимо только к квадратным матрицам, то SVD позволяет разложить абсолютно любую прямоугольную матрицу. Эта математическая процедура стала одним из важнейших открытий XX века, заложив фундамент для современной науки о данных, обработки изображений и систем искусственного интеллекта.

Геометрический смысл сингулярного разложения потрясающе красив. Любое линейное преобразование (которое описывается матрицей) можно разбить на три простых последовательных шага. Первый шаг — это поворот исходного пространства. Второй шаг — масштабирование (растяжение или сжатие) пространства вдоль новых координатных осей. Третий шаг — еще один поворот. Матрица SVD буквально предоставляет нам матрицы этих поворотов (ортогональные матрицы сингулярных векторов) и коэффициенты растяжения (диагональную матрицу сингулярных чисел). Чем больше сингулярное число, тем больше информации или «энергии» несет соответствующее направление.

Сжатие информации и метод главных компонент (PCA)

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

Этот принцип лежит в основе Метод Главных Компонент (PCA) — важнейшего алгоритма снижения размерности данных. Например, черно-белое изображение лица размером 1000x1000 пикселей (миллион чисел) с помощью усеченного SVD можно сжать до передачи всего 50 векторов, и глаз человека не заметит разницы в качестве. В знаменитом конкурсе Netflix Prize с призовым фондом в миллион долларов, победу одержала команда, алгоритм которой базировался именно на SVD. Матрица, где строками были пользователи, а столбцами — фильмы, содержала огромное количество пропусков (никто не смотрел все фильмы). SVD позволило выявить скрытые (латентные) факторы вкусов пользователей и жанров фильмов, блестяще предсказывая оценки, которые пользователи поставили бы еще не просмотренным картинам.

Подробнее

Метод конечных объемов (МКО): фундамент вычислительной гидродинамики

Законы сохранения как основа физики сплошных сред

Ранее мы рассматривали метод конечных разностей (МКР) и метод конечных элементов (МКЭ). Однако, когда дело доходит до моделирования течений газов и жидкостей (аэродинамика самолетов, течения в газопроводах, взрывы, горение в двигателях), инженеры обращаются к третьему, специализированному инструменту — методу конечных объемов (МКО, или Finite Volume Method, FVM). Именно этот алгоритм заложен в основу большинства ведущих коммерческих и открытых пакетов вычислительной гидродинамики (CFD), таких как ANSYS Fluent, CFX и OpenFOAM.

Причина такого выбора кроется в самой сути дифференциальных уравнений газовой динамики (системы уравнений Эйлера или Навье-Стокса). Эти уравнения являются математической записью фундаментальных законов природы: сохранения массы, сохранения импульса и сохранения энергии. МКР и МКЭ оперируют узлами сетки и пытаются аппроксимировать дифференциальные операторы. При грубой сетке они могут математически «потерять» часть массы или энергии. Метод конечных объемов, напротив, строится не на дифференциальных, а на интегральных формах законов сохранения, гарантируя строгий баланс физических величин на любой, даже самой грубой сетке.

Контрольные объемы и расчет потоков на гранях

Суть МКО состоит в следующем: расчетная область со сложной геометрией разбивается на множество мелких, непересекающихся ячеек (многогранников). Эти ячейки называются контрольными объемами. Искомые величины (плотность, давление, скорость) хранятся не в узлах, а в центрах этих объемов и представляют собой средние значения внутри ячейки.

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

Улавливание ударных волн и решатели Римана

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

Для решения этой проблемы российский математик Сергей Годунов в 1959 году предложил гениальную идею. Он предложил рассматривать разрыв параметров на гране между двумя ячейками как локальную задачу о распаде произвольного разрыва (задачу Римана в газовой динамике). Решая задачу Римана на каждой грани, алгоритм получает точные, физически обоснованные потоки, которые естественным образом учитывают направление распространения волн (скорость звука). Развитие этой идеи привело к созданию современных высокоразрешающих схем с ограничителями потоков (TVD-схемы, алгоритмы Roe, HLLC). Благодаря им МКО способен рассчитывать сложнейшие сверхзвуковые течения с резкими фронтами ударных волн без численных осцилляций, обеспечивая человечество инструментами для проектирования гиперзвуковых аппаратов и турбин.

Подробнее

Оптимизация в машинном обучении: стохастический градиентный спуск (SGD)

Большие данные и паралич классических методов

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

Проблема заключается в масштабах. Функция потерь нейросети складывается из ошибок на каждом отдельном примере из обучающей выборки. Если ваша выборка состоит из 10 миллионов изображений (как в датасете ImageNet), то для того, чтобы сделать всего один крошечный шаг градиентного спуска, компьютеру нужно прогнать через огромную нейросеть все 10 миллионов картинок, вычислить ошибку для каждой, сложить их и найти производные. Это чудовищно долго. Классический (полнопакетный) градиентный спуск в эпоху Big Data оказался вычислительно парализован.

Стохастичность: скорость в обмен на шум

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

Это кардинально меняет картину. Вместо того чтобы ждать часы ради одного абсолютно точного шага, нейросеть делает тысячи приблизительных шагов в секунду. Да, каждый отдельный шаг SGD сильно зашумлен — градиент от одного примера может указывать не совсем в сторону глобального минимума, заставляя траекторию поиска сильно «дрожать» и петлять. Однако, согласно закону больших чисел, математическое ожидание (среднее направление) этих хаотичных скачков верно ведет алгоритм к минимуму функции. Более того, доказано, что этот стохастический шум работает как форма регуляризации: он помогает алгоритму «выскакивать» из мелких локальных минимумов, что крайне полезно в невыпуклом ландшафте функций потерь глубоких сетей.

Мини-батчи и адаптивные оптимизаторы (Adam)

В чистом виде SGD (по одному примеру) плохо использует возможности современных видеокарт (GPU), которые созданы для параллельных матричных вычислений. Поэтому золотым стандартом индустрии стал компромисс — SGD по мини-батчам (Mini-batch SGD). Выборка делится на небольшие случайные пакеты (от 32 до 512 примеров). Градиент вычисляется по этому пакету. Это дает и достаточную скорость благодаря матричным вычислениям на GPU, и спасительный стохастический шум.

Чтобы справиться с проблемой подбора скорости обучения (learning rate) и зигзагообразными колебаниями SGD в оврагах, математики разработали мощные надстройки. Введение «импульса» (Momentum) позволяет алгоритму накапливать скорость, двигаясь в постоянном направлении, и пробивать мелкие препятствия по инерции. Вершиной этой эволюции стал алгоритм Adam (Adaptive Moment Estimation), который для каждого отдельного из миллионов параметров нейросети вычисляет свою собственную, адаптивную скорость обучения, базируясь на скользящих средних градиентов и их квадратов. Сегодня Adam и его модификации (AdamW) являются инструментами по умолчанию для обучения подавляющего большинства систем искусственного интеллекта, от распознавания образов до больших языковых моделей.

Подробнее

Оценивание состояния систем: фильтр Калмана и анализ зашумленных данных

Проблема скрытых состояний и неточных датчиков

В задачах управления сложными динамическими системами — будь то беспилотный автомобиль, межконтинентальная ракета или экономическая модель — мы сталкиваемся с двумя фундаментальными проблемами. Во-первых, мы не всегда можем измерить желаемые параметры напрямую (скрытое состояние). Например, мы можем измерять температуру снаружи двигателя, но нам нужно знать температуру в его центре. Во-вторых, любые измерения в реальном мире подвержены случайным ошибкам и шумам датчиков. Данные GPS скачут на метры, акселерометры накапливают ошибку дрейфа. Как из этого хаоса неточных, зашумленных данных извлечь истинное состояние системы в режиме реального времени?

Решение этой задачи было найдено в 1960 году Рудольфом Калманом. Его алгоритм, получивший название Фильтр Калмана, стал одним из важнейших открытий в истории вычислительной математики и теории управления. Именно этот фильтр позволил бортовым компьютерам миссии «Аполлон» точно оценивать траекторию спускаемого модуля и успешно приземлиться на Луну, несмотря на сильнейшие шумы радаров.

Математика предсказаний и коррекций

Фильтр Калмана работает по принципу рекурсивного байесовского оценивания. Он не требует хранения огромного массива прошлых данных; ему нужны только результаты предыдущего шага и текущее измерение. Алгоритм элегантно объединяет два источника информации: физическую (математическую) модель системы и показания датчиков. Работа фильтра состоит из бесконечно повторяющегося двухтактного цикла: Предсказание (Predict) и Коррекция (Update).

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

Магия коэффициента усиления Калмана

Суть алгоритма скрыта в вычислении специальной матрицы — коэффициента усиления Калмана (Kalman Gain). Этот коэффициент определяет, кому алгоритм должен доверять больше в данный момент: своей математической модели или показаниям датчика. Если модель очень точная, а датчик сильно шумит (например, во время бури), коэффициент усиления уменьшается, и фильтр игнорирует скачки показаний, опираясь на физику. Если же модель неточна (возник непредсказуемый порыв ветра), а датчик высокоточный, коэффициент усиления возрастает, и фильтр жестко корректирует свою внутреннюю оценку по данным прибора.

Путем умножения матриц и обновления ковариаций, фильтр Калмана математически доказывает, что его оценка является оптимальной (минимизирует среднеквадратичную ошибку) для линейных систем с нормальным гауссовским шумом. Для нелинейных систем (с которыми мы работаем в 99% случаев) был разработан Расширенный фильтр Калмана (EKF), который на каждом шаге линеаризует уравнения с помощью матриц Якоби, оставаясь сегодня абсолютным стандартом в робототехнике, навигации дронов и системах автопилотов.

Подробнее

Глобальная оптимизация: генетические алгоритмы и эволюционные вычисления

Ловушка локальных минимумов

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

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

Биологическая метафора: хромосомы, популяции и отбор

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

На каждом шаге (в каждом поколении) алгоритм оценивает приспособленность (fitness) каждой особи. Приспособленность — это просто значение целевой функции, которую мы хотим минимизировать или максимизировать. Чем лучше решение, тем выше шансы особи на выживание. Затем вступает в дело оператор селекции (рулетка или турнирный отбор): алгоритм с большей вероятностью выбирает лучших особей для формирования «родительского пула», безжалостно отсеивая слабые, неудачные решения. Это математическое воплощение принципа выживания наиболее приспособленных.

Кроссовер и мутация: двигатели математического прогресса

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

Однако селекция и кроссовер могут привести к тому, что популяция быстро выродится и станет состоять из одинаковых особей, застряв в локальном экстремуме. Чтобы поддерживать генетическое разнообразие и исследовать новые, неизведанные участки пространства, применяется оператор мутации. С небольшой вероятностью (например, 1%) алгоритм случайно изменяет один или несколько генов (битов) у потомков. Мутация позволяет алгоритму случайно «перепрыгивать» через холмы целевой функции, вырываясь из ловушек локальных минимумов. Повторяя эти циклы (оценка, селекция, скрещивание, мутация) на протяжении сотен поколений, генетический алгоритм стягивает популяцию к глобальному оптимуму даже для недифференцируемых, разрывных и зашумленных функций, с которыми не справляется ни один градиентный метод.

Подробнее

Численное решение краевых задач: метод стрельбы и метод прогонки

Отличие краевой задачи от задачи Коши

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

Краевые задачи значительно сложнее вычислительно. Мы не можем просто взять метод Рунге-Кутты и начать интегрировать от одного конца к другому, потому что нам не хватает начальных данных (мы знаем координату, но не знаем производную). Приходится применять специальные алгоритмы, которые ищут решение, удовлетворяющее одновременно всем граничным условиям на обоих концах интервала.

Метод стрельбы: искусство прицеливания в математике

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

Алгоритм действует так: мы выбираем некоторый случайный или приближенный начальный угол (значение производной) и решаем систему как обычную задачу Коши методами Эйлера или Рунге-Кутты. Дойдя до конца интервала, мы смотрим, куда попал наш «снаряд». Скорее всего, мы промахнулись мимо требуемого граничного условия. Мы вычисляем величину промаха (невязку). Затем мы меняем начальный угол, стреляем еще раз и снова оцениваем промах. Задача сводится к поиску корня нелинейного уравнения: нам нужно подобрать такой начальный параметр, при котором невязка на правом конце станет равной нулю. Для быстрого поиска этого параметра обычно применяется метод Ньютона или метод секущих.

Метод конечных разностей и алгоритм прогонки

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

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

Подробнее

Дискретное вейвлет-преобразование: альтернатива Фурье в анализе сигналов

Анализ сигналов: выход за пределы синусоид

На протяжении почти двух столетий непревзойденным королем спектрального анализа оставалось преобразование Фурье. Его суть заключается в разложении любого сложного сигнала на сумму бесконечных, периодических, идеально гладких синусоид и косинусоид разных частот и амплитуд. Алгоритм Быстрого Преобразования Фурье (БПФ) сделал цифровую обработку звука и изображений реальностью. Однако классическое преобразование Фурье имеет один фундаментальный недостаток: оно глобально. Синусоида длится от минус бесконечности до плюс бесконечности. Фурье-анализ идеально показывает, какие частоты присутствуют в сигнале, но он абсолютно не способен сказать, в какой именно момент времени эта частота появилась или исчезла.

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

Вейвлеты: маленькие волны с конечной энергией

Этим новым инструментом стало вейвлет-преобразование, активно развивавшееся в 1980-х годах (труды Морле, Гроссмана, Ив Мейер и Ингрид Добеши). Слово «вейвлет» (wavelet) переводится как «маленькая волна» или «всплеск». В отличие от синусоиды, вейвлет отличен от нуля только на коротком отрезке времени, а за его пределами быстро затухает до нуля.

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

Дискретное преобразование и индустриальные стандарты

На практике в вычислительной технике используется Дискретное Вейвлет-Преобразование (ДВП). В нем масштабы и сдвиги изменяются не плавно, а по степеням двойки (диадическая сетка). Алгоритм Малла (Stephane Mallat) позволил вычислять ДВП с помощью каскада цифровых фильтров: на каждом шаге сигнал пропускается через низкочастотный и высокочастотный фильтры, разделяя данные на аппроксимирующие (грубые) и детализирующие коэффициенты.

Благодаря своей способности концентрировать энергию сигнала в небольшом количестве коэффициентов (большинство коэффициентов ДВП для реальных сигналов оказываются близки к нулю), вейвлеты совершили революцию в сжатии данных. Именно вейвлет-преобразование (конкретно, вейвлеты Коэна-Добеши-Фово) лежит в основе графического стандарта JPEG 2000, обеспечивая сжатие изображений высокого разрешения без характерных для обычного JPEG раздражающих квадратных «артефактов блочности». Также вейвлеты стали золотым стандартом для очистки медицинских и акустических сигналов от шумов (алгоритмы wavelet denoising/shrinkage).

Подробнее

Итерационные методы подпространств Крылова для гигантских матриц

Как решать системы с миллиардами уравнений?

При дискретизации трехмерных дифференциальных уравнений в частных производных (например, в аэродинамике или гидрогеологии) методом конечных элементов возникают колоссальные системы линейных алгебраических уравнений (СЛАУ). Размерность таких матриц может достигать сотен миллионов и даже миллиардов строк. Использование прямого метода Гаусса для таких гигантов невозможно: даже суперкомпьютеру не хватит оперативной памяти для хранения матрицы, не говоря уже о годах, необходимых для расчетов.

К счастью, такие матрицы являются крайне разреженными. Более 99% их элементов — это нули. Для их хранения используются специальные форматы (например, CSR — Compressed Sparse Row), в которых хранятся только ненулевые элементы и их координаты. Однако прямые методы в процессе исключения неизбежно превращают нули в ненулевые элементы (проблема заполнения). Поэтому для решения таких огромных разреженных СЛАУ применяются исключительно современные итерационные методы, вершиной эволюции которых являются алгоритмы на основе подпространств Крылова.

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

В 1931 году русский математик Алексей Крылов опубликовал работу, которая заложила основы нового класса методов. Идея состоит в следующем: поскольку мы можем быстро умножать нашу гигантскую разреженную матрицу (назовем ее A) на вектор, почему бы не использовать операцию умножения для построения базиса, в котором мы будем искать решение?

Если мы возьмем начальный вектор невязки (r) и начнем последовательно умножать его на матрицу A, мы получим последовательность векторов: r, Ar, A^2r, A^3r и так далее. Линейная оболочка (пространство, натянутое на эти векторы) называется подпространством Крылова. Гениальность подхода заключается в том, что вместо решения системы в исходном пространстве с миллионами измерений, мы ищем приближенное решение (проекцию) в подпространстве Крылова, размерность которого мала (например, 50 или 100). На каждом шаге (итерации Арнольди или Ланцоша) к подпространству добавляется новый ортогональный вектор, и приближение становится все точнее.

GMRES и Метод сопряженных градиентов (CG)

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

Однако если матрица несимметричная (что часто бывает в задачах конвекции-диффузии и гидродинамики), метод CG не работает. Для таких матриц в 1986 году Саадом и Шульцем был разработан метод GMRES (Generalized Minimal Residual method). Метод GMRES минимизирует норму невязки на всем подпространстве Крылова. В отличие от CG, он вынужден хранить в памяти все базисные векторы, поэтому на практике применяется с «рестартами» (GMRES(m)): после достижения m итераций накопленные векторы сбрасываются, и процесс начинается заново с текущего наилучшего приближения. Сегодня алгоритмы подпространств Крылова, оснащенные мощными предобуславливателями (preconditioners) типа ILU или Algebraic Multigrid (AMG), являются абсолютным индустриальным стандартом при суперкомпьютерном моделировании.

Подробнее

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

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

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

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

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

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

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

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

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

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

Подробнее
Subscribe to this RSS feed

ГДЗ по математике

1 класс 2 класс 3 класс
4 класс 5 класс 6 класс

Книги

Фильмы

Старый сайт

Соц. сети