Main menu
Численные методы

Численные методы (100)

Предобуславливание (Preconditioning): ключ к решению гигантских СЛАУ

Узкое место итерационных методов Крылова

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

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

Математическая трансформация системы

Для преодоления этого вычислительного тупика применяется предобуславливание (Preconditioning). Идея заключается в преобразовании исходной СЛАУ (Ax = b) в другую, математически эквивалентную систему, которая имеет точно такое же решение, но гораздо лучшее число обусловленности. Для этого вводится матрица предобуславливателя (обозначим ее M).

Мы умножаем левую и правую части системы на обратную матрицу предобуславливателя: M^(-1)Ax = M^(-1)b. Теперь алгоритм Крылова решает новую систему с матрицей M^(-1)A. Идеальным предобуславливателем была бы сама матрица A (тогда M^(-1)A превратилась бы в единичную матрицу, и алгоритм сошелся бы за 1 итерацию). Но обращение A — это именно то, чего мы избегаем. Следовательно, матрица M должна удовлетворять двум противоречивым требованиям: она должна быть «максимально похожа» на матрицу A (чтобы сгруппировать собственные значения), и при этом решение системы My = z должно вычисляться алгоритмически очень быстро и дешево (так как это придется делать на каждой итерации).

Популярные алгоритмы: ILU и диагональное масштабирование

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

Золотым стандартом в индустрии стало Неполное LU-разложение (Incomplete LU factorization, ILU). Алгоритм пытается сделать стандартное разложение матрицы A на нижнюю и верхнюю треугольные матрицы, но принудительно отбрасывает (заменяет нулями) те новые элементы (fill-ins), которые возникают в процессе исключения. В результате мы получаем приближенные матрицы L и U, которые сохраняют идеальную разреженность исходной системы. Решение системы с такими треугольными матрицами выполняется мгновенно прямой и обратной подстановкой. Баланс между точностью ILU (задаваемой порогом отбрасывания) и скоростью одной итерации является главным искусством настройки современных коммерческих решателей. Для самых тяжелых задач применяются многосеточные предобуславливатели (AMG), которые обеспечивают теоретическую линейную масштабируемость при решении миллиардных систем.

Подробнее

Метод сглаженных частиц (SPH): бессеточная вычислительная гидродинамика

Проблема подвижных границ и сильных деформаций в сеточных методах

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

Попытка отслеживать свободную поверхность (поверхность раздела газ-жидкость) на неподвижной сетке требует использования сложных алгоритмов типа Volume of Fluid (VOF) или Level Set, которые размазывают границу и подвержены численной диффузии. Использование же деформируемых (лагранжевых) сеток приводит к их быстрому скручиванию и вырождению ячеек, требуя постоянного, вычислительно дорогого перестроения (ремешинга). Инженеры и физики остро нуждались в методе, который вообще не зависит от геометрической сетки.

Лагранжев подход: жидкость как ансамбль частиц

В 1977 году Люси, а также Гинголд и Монаган, независимо друг от друга предложили метод сглаженных частиц (Smoothed Particle Hydrodynamics, SPH) для решения задач астрофизики. SPH — это чисто лагранжев бессеточный метод. Среда (жидкость, газ или деформируемое твердое тело) представляется в виде набора дискретных макроскопических частиц. Каждая частица обладает фиксированной массой и несет в себе информацию о плотности, давлении, внутренней энергии и скорости.

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

Ядро сглаживания и аппроксимация производных

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

Весовые коэффициенты определяются специальной колоколообразной функцией — ядром сглаживания (Smoothing kernel). Ядро имеет радиус действия (smoothing length). Частицы, находящиеся за пределами этого радиуса, не влияют на расчет, что позволяет использовать алгоритмы быстрого поиска соседей (например, KD-деревья). Самое важное свойство SPH заключается в том, что пространственная производная от физической величины аналитически переносится на саму функцию ядра! То есть, чтобы найти градиент давления для уравнения импульса, компьютеру не нужно брать разности между частицами; он просто берет аналитическую производную от известной функции ядра. Благодаря своей гибкости метод SPH сегодня является стандартом не только в моделировании баллистики и цунами, но и в компьютерной графике для создания реалистичной воды в голливудских фильмах.

Подробнее

Алгоритмы минимизации с ограничениями: метод штрафных функций и барьерные методы

Реальный мир диктует жесткие условия

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

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

Метод штрафных функций (Exterior Penalty Method)

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

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

Барьерные методы (Interior Point Methods)

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

Эту проблему решают барьерные методы (методы внутренних точек). Они действуют зеркально: алгоритм всегда строго удерживается внутри разрешенной зоны. К целевой функции прибавляется барьерная функция (часто это логарифм от расстояния до границы). В центре разрешенной области барьер не мешает поиску. Но как только алгоритм пытается приблизиться к опасной границе, барьерная функция устремляется в бесконечность (строит математическую бетонную стену), не позволяя нарушить физические законы. Постепенно ослабляя вес барьера (уменьшая барьерный параметр), алгоритм позволяет решению медленно и безопасно «доползти» до оптимальной точки на границе. Современные алгоритмы прямо-двойственных внутренних точек (Primal-Dual Interior Point Methods), развившие эту концепцию, произвели фурор в 1990-х годах, став абсолютным стандартом решения гигантских задач линейного и нелинейного программирования в логистике, энергетике и финансах.

Подробнее

Псевдоспектральные методы решения нелинейных эволюционных уравнений

Нелинейные волны: от солитонов до оптоволокна

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

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

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

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

Гениальным выходом из ситуации стали псевдоспектральные методы (Pseudospectral methods, или методы Фурье-коллокации). Их философия проста: «Выполняй каждую математическую операцию в том пространстве, где она алгоритмически дешевле всего». В этом методе вычисление пространственных производных происходит в частотном пространстве (где это тривиальное и точное умножение на волновое число). Но как только дело доходит до нелинейного члена, алгоритм вызывает Быстрое преобразование Фурье (БПФ) и мгновенно переносит данные в физическое пространство!

Расщепление по физическим процессам (Split-Step)

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

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

Подробнее

Численные методы в задачах фильтрации: закон Дарси и уравнение Маскета

Движение флюидов в пористых средах

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

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

Многофазная фильтрация и модель Бакли-Леверетта

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

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

Метод IMPES: разделяй и властвуй в гидродинамике пласта

Золотым стандартом в коммерческих симуляторах резервуаров (таких как Eclipse или tNavigator) является численный алгоритм IMPES (Implicit Pressure Explicit Saturation — неявное давление, явная насыщенность). Идея алгоритма основана на физическом факте: давление в горных породах (звуковая волна) распространяется в тысячи раз быстрее, чем физически перемещаются фронты жидкостей (насыщенность).

Алгоритм IMPES искусственно разделяет (расщепляет) сложную систему на две независимые части. На первом микрошаге алгоритм считает насыщенности всех фаз «замороженными» (константами). Он составляет и неявно решает гигантскую СЛАУ для нахождения нового поля давлений. Так как схема неявная, она абсолютно устойчива и позволяет брать большие шаги по времени. Получив точное распределение давления (а значит, и скоростей потоков), на втором этапе алгоритм по явной формуле вычисляет перенос масс жидкостей (новые насыщенности) вдоль этих потоков. Для устойчивого расчета насыщенности критически важно использовать методы направленных разностей (Upwind schemes) с ограничителями потоков, чтобы избежать нефизичных осцилляций на фронтах вытеснения. Численное моделирование пластов (Reservoir Simulation) — это одна из тех областей, где вычислительная математика напрямую генерирует миллиардные прибыли для корпораций.

Подробнее

Метод фиктивных областей в задачах со сложной геометрией

Сложность построения сеток для криволинейных границ

В вычислительной математике, особенно при решении краевых задач для дифференциальных уравнений в частных производных (УЧП), геометрия физического объекта диктует метод решения. Если мы моделируем теплопроводность в идеальном металлическом кубе, достаточно использовать тривиальную прямоугольную (декартову) сетку. Матрицы, получаемые на таких сетках, обладают идеальной структурой, а системы уравнений решаются мгновенно с помощью быстрых алгоритмов (например, многосеточных методов или БПФ).

Но что делать, если геометрия детали невероятно сложна, содержит внутренние полости, острые углы и криволинейные поверхности? Классический путь — это использование неструктурированных треугольных или тетраэдральных сеток (в методе конечных элементов). Однако процесс генерации качественной 3D-сетки вокруг сложной детали (например, блока цилиндров двигателя) является тяжелейшей задачей, которая может занимать дни работы инженера и часто требует ручного вмешательства для устранения «вырожденных» элементов. Метод фиктивных областей (Fictitious Domain Method) предлагает радикальную и математически красивую альтернативу: вместо того чтобы подстраивать сетку под сложную форму тела, мы подстроим сами уравнения.

Погружение сложной детали в простую «коробку»

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

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

Метод штрафа: математическое моделирование жесткой стены

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

Наличие такого гигантского множителя заставляет уравнение в фиктивной области «вырождаться». Чтобы математический баланс не нарушился, алгоритм вынужден жестко приравнять искомую температуру в фиктивной области к нулю (или другому заданному краевому значению). Таким образом, мы получаем математическую модель, в которой фиктивная область действует как абсолютно жесткая «ледяная глыба», заставляющая решение подчиняться краевым условиям прямо на криволинейной границе детали (которая проходит сквозь прямоугольную сетку). Метод фиктивных областей позволяет использовать сверхбыстрые регулярные решатели для задач любой формы, став основой для многих современных CAD-пакетов, не требующих от инженера ручного построения сеток (Meshless-подобные подходы).

Подробнее

Численное решение уравнений в свертках и алгоритмы деконволюции

Искажение сигналов и размытие изображений

В физике, астрономии, оптике и цифровой обработке сигналов мы постоянно сталкиваемся с тем, что измерительная аппаратура несовершенна. Любой прибор (будь то телескоп, микроскоп или сейсмограф) неизбежно искажает истинный физический сигнал. Математически процесс искажения линейной стационарной системой (LTI-системой) описывается операцией свертки: истинный сигнал «сворачивается» с так называемой аппаратной функцией прибора (функцией рассеяния точки, PSF). В результате мы получаем размытую, сглаженную картину реальности.

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

Деление в частотной области и фильтр Винера

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

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

Итерационная деконволюция Ричардсона-Люси

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

Алгоритм Ричардсона-Люси работает итеративно: на каждом шаге он берет текущее приближение истинного изображения, сворачивает его с аппаратной функцией (симулируя искажение), сравнивает полученный результат с реальной размытой фотографией путем деления, а затем использует эту разность для коррекции приближения на следующем шаге. Гениальность метода состоит в том, что он математически гарантирует неотрицательность пикселей на любой итерации и отлично сохраняет полную энергию сигнала. Именно алгоритм Ричардсона-Люси стал спасением для космического телескопа «Хаббл» в 1990 году, когда выяснилось, что его главное зеркало отполировано с ошибкой (имеет сферическую аберрацию): численные методы деконволюции позволили ученым получать резкие снимки галактик еще до прибытия ремонтной экспедиции с корректирующей оптикой.

Подробнее

Методы фиктивных областей и погруженных границ (Immersed Boundary Method)

Проклятие подвижных сеток в гидродинамике

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

Если использовать классический метод конечных объемов, нам придется на каждом крошечном шаге по времени заново перестраивать всю вычислительную сетку (ремешинг), чтобы она подстраивалась под новые границы объекта. Перестроение объемной 3D-сетки — это невероятно медленная, тяжелая алгоритмическая задача, которая съедает до 90% всего процессорного времени, а сильные деформации ячеек быстро разрушают точность вычислений. Для решения этой фундаментальной проблемы механики сплошных сред был разработан радикальный подход — отказ от согласования сеток. Этот подход известен как Метод фиктивных областей (Fictitious Domain Method) и его знаменитая вариация — Метод погруженных границ (Immersed Boundary Method, IBM).

Идея метода погруженных границ Чарльза Пескина

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

Первая сетка — эйлерова. Она представляет собой простую, неподвижную, равномерную и строгую прямоугольную (картезианскую) сетку, натянутую на весь объем пространства. В ней решаются уравнения течения жидкости (Навье-Стокса). Она никогда не меняется, что позволяет использовать сверхбыстрые решатели на базе БПФ. Вторая сетка — лагранжева. Это эластичная, подвижная сетка, состоящая только из одномерных кривых или двумерных поверхностей, которая описывает само движущееся тело (например, створку клапана). Эта сетка «плавает» сквозь неподвижную жидкостную сетку, совершенно не обращая внимания на ее узлы.

Взаимодействие через дельта-функции Дирака

Главный математический вызов метода — как заставить неподвижную жидкость почувствовать, что сквозь нее летит твердое тело или мембрана? Здесь на помощь приходит мощный инструмент функционального анализа — дельта-функция Дирака.

Поверхностная лагранжева сетка (тело) вычисляет упругие силы, возникающие при ее деформации. Затем эти поверхностные силы с помощью специальной сглаженной дельта-функции «размазываются» (интерполируются) на ближайшие узлы неподвижной жидкостной сетки, превращаясь в объемные источники силы в уравнениях Навье-Стокса. Жидкость, получив этот силовой толчок, меняет свое поле скоростей. Затем жидкость диктует свою новую скорость точкам мембраны, заставляя их сместиться на новый шаг по времени. Этот изящный алгоритм взаимодействия (Fluid-Structure Interaction, FSI) исключает необходимость перестроения сеток. Сегодня методы погруженных границ являются мощнейшим инструментом для моделирования гемодинамики, движения клеточных мембран, создания цифровых двойников кровеносной системы человека и расчетов сложной многофазной промышленной химии.

Подробнее

Метод сплайн-коллокации: гладкие решения краевых задач без интегрирования

Альтернатива Галеркину: сила точечного совпадения

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

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

Базис B-сплайнов: локальность и высокая гладкость

Теоретически в качестве пробной функции можно взять один длинный глобальный полином (многочлен). Но, как мы помним из феномена Рунге, глобальные многочлены высоких степеней сильно осциллируют и ведут себя нестабильно. Идеальным базисом для метода коллокации оказались B-сплайны (базисные сплайны), разработанные Исааком Шёнбергом и Карлом де Бором.

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

Оптимальные узлы: корни полиномов Чебышева и Гаусса

Самым критическим моментом в методе сплайн-коллокации является выбор расположения самих точек коллокации. Если расставить эти точки равномерно, метод может потерять устойчивость или сходиться очень медленно. Математический анализ, проведенный Карлом де Бором и Ричардом Сварцем в 1973 году, доказал потрясающую вещь: точность метода кардинально возрастает (явление суперсходимости), если в качестве точек коллокации внутри каждого сеточного интервала использовать корни ортогональных полиномов (узлы Гаусса-Лежандра или Чебышева).

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

Подробнее

Численное интегрирование функций с сингулярностями: метод Tanh-Sinh

Когда подынтегральная функция уходит в бесконечность

Классические методы численного интегрирования (квадратуры Ньютона-Котеса, формулы Симпсона или Гаусса) блестяще справляются с гладкими функциями, которые заданы на конечном отрезке. Однако в физических задачах математикам постоянно приходится сталкиваться с интегралами, которые имеют математические особенности — сингулярности. Например, при расчетах потенциала заряженной нити или в задачах механики разрушения (напряжения на кончике трещины) подынтегральная функция может уходить в бесконечность на краях отрезка или иметь вертикальные асимптоты. Классический алгоритм, столкнувшись с бесконечностью (делением на ноль), выдаст ошибку Not-a-Number (NaN) или будет сходиться катастрофически медленно.

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

Магия нелинейных замен переменных

Второй, гораздо более мощный алгоритмический подход основан на нелинейных заменах координат. Идея заключается в том, чтобы так исказить ось интегрирования, чтобы «растянуть» проблемные участки (около сингулярности) и свести бесконечности к нулю. Одной из первых успешных замен такого рода стала подстановка Игараси-Мори-Такахаси (IMT).

Но настоящим прорывом и шедевром японской вычислительной школы стало создание в 1974 году квадратурной формулы гиперболического тангенса-синуса (Tanh-Sinh quadrature, или метод двойной экспоненциальной интеграции). Этот метод был разработан Хидетаси Такахаси и Масатаке Мори. Идея заключается во введении сложнейшей подстановки: x = tanh(pi/2 * sinh(t)). Эта функция обладает потрясающим топологическим свойством: она отображает бесконечный интервал (от минус до плюс бесконечности) строго в конечный отрезок [-1, 1].

Двойное экспоненциальное затухание

При выполнении этой замены дифференциал интегрирования (якобиан) порождает функцию, которая убывает к нулю с феноменальной скоростью — по закону двойной экспоненты: exp(-C * exp(|t|)). Эта скорость убывания настолько сокрушительна, что она способна мгновенно подавить (помножить на ноль) абсолютно любую математическую сингулярность исходной функции на краях отрезка!

После такой замены переменной мы получаем новый, гладкий и колоколообразный интеграл на бесконечной прямой, для вычисления которого применяется самое простое правило трапеций (да, базовая формула трапеций на бесконечном интервале для гладких убывающих функций работает точнее сложных методов Гаусса). Алгоритм Tanh-Sinh позволяет с феноменальной легкостью брать интегралы с любыми концевыми сингулярностями, логарифмическими разрывами или бесконечными производными, достигая невероятной точности в тысячи знаков после запятой в системах произвольной разрядности. Сегодня этот алгоритм встроен в ядра таких мощных математических систем, как MPMath в Python и Boost.Math в C++.

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

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

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

Книги

Фильмы

Старый сайт

Соц. сети