Высшая математика

Решение уравнений Навье-Стокса: алгоритмы SIMPLE и PISO для несжимаемых жидкостей

Вычислительный парадокс давления в несжимаемых потоках

Математическое моделирование турбулентных течений воздуха вокруг фюзеляжей самолетов или потоков вязкой воды в промышленных трубах базируется на интегрировании сложнейших уравнений Навье-Стокса. Для подавляющего большинства земных инженерных задач обычные жидкости и газы при относительно низких скоростях (до числа Маха, равного 0.3) считаются физически несжимаемыми: их термодинамическая плотность строго постоянна. Но именно эта кажущаяся физическая простота порождает колоссальную математическую проблему при численном решении системы на компьютере. Уравнения импульса (закон Ньютона) легко позволяют нам найти векторное поле скоростей потока, но только если нам заранее известно скалярное поле внутреннего давления. Но как математически найти само давление?

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

Связка давления и скорости: алгоритм предиктор-корректор SIMPLE

Для блестящего решения этой нетривиальной проблемы связи давления и скорости (pressure-velocity coupling) в 1972 году физики Сполдинг и Патанкар разработали легендарный алгоритм SIMPLE (Semi-Implicit Method for Pressure Linked Equations). Идея этого численного алгоритма концептуально основана на классическом методе предиктор-корректор.

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

Подробнее

Методы доверительных областей (Trust Region) в нелинейной оптимизации

Мощная альтернатива классическому линейному поиску

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

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

Построение суррогатной квадратичной модели и адаптивный радиус

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

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

Подробнее

Численная оптимизация топологии: как математика создает бионический дизайн

От интуиции инженера к строгой вычислительной математике

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

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

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

Самым распространенным и вычислительно эффективным алгоритмом топологической оптимизации в мировой индустрии является метод SIMP (Solid Isotropic Material with Penalization). Алгоритм начинает работу с того, что разбивает всю доступную расчетную область на мелкие кубические конечные элементы (воксели). Изначально каждому такому элементу программно присваивается математическая «плотность» от 0 (полная воздушная пустота) до 1 (цельный монолитный металл). Проблема заключается в том, что в реальном физическом мире материал детали не может быть полупрозрачным, с плотностью 0.5. Нам на выходе нужен строгий дискретный ответ: либо там есть сталь, либо там воздух.

Чтобы заставить математический алгоритм стремиться к четким черным и белым границам, метод SIMP вводит штрафование (Penalization) промежуточных серых плотностей. Жесткость каждого элемента привязывается к его расчетной плотности через нелинейную степенную функцию (обычно плотность возводится в куб). Таким образом, полупустой элемент со значением 0.5 будет иметь всего 12.5% от максимальной жесткости, но «весить» как половина полного. Оптимизатору (обычно используется метод MMA - Method of Moving Asymptotes) становится математически крайне невыгодно и энергетически дорого оставлять такие серые зоны. Он агрессивно выталкивает плотность каждого элемента к строгим краям: либо в 0, либо в 1. На каждом шаге итерации алгоритм решает гигантскую СЛАУ упругости, вычисляет градиенты функции чувствительности и плавно обновляет плотности элементов. С развитием промышленной 3D-печати этот метод стал абсолютным стандартом в аэрокосмической отрасли.

Подробнее

Метод граничных элементов (МГЭ): снижение размерности задачи на единицу

Снижение размерности: от объема к поверхности

Метод конечных элементов (МКЭ) и метод конечных разностей (МКР) безусловно доминируют в современных инженерных расчетах. Однако оба эти метода являются полностью объемными: они требуют построения густой расчетной сетки во всей 3D-области математического моделирования. Если вам нужно рассчитать обтекание подводной лодки в открытом океане или напряжение вокруг микротрещины в огромном стальном блоке, объемная сетка съест колоссальное количество памяти суперкомпьютера. В таких специфических случаях на помощь приходит элегантный математический аппарат — Метод граничных элементов (МГЭ, или Boundary Element Method, BEM).

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

Интегральные уравнения и фундаментальные решения

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

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

Подробнее

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

От классических полиномов к глубокому обучению

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

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

Физически информированные нейросети (PINNs)

Настоящей революцией в численном решении дифференциальных уравнений в частных производных (УЧП) стало появление физически информированных нейронных сетей (Physics-Informed Neural Networks, PINNs). Идея PINN заключается во внедрении физических законов (самих дифференциальных уравнений) непосредственно в функцию потерь (Loss function) нейронной сети во время ее математического обучения.

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

Подробнее

Разреженные матрицы и их хранение: форматы CSR, CSC и COO

Экономия памяти в гигантских системах

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

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

Координатный формат (COO)

Самым интуитивно понятным и простым способом сохранения разреженной матрицы является координатный формат (Coordinate format, COO). В этом формате матрица хранится в виде трех одномерных массивов одинаковой длины (равной количеству ненулевых элементов, NNZ). Первый массив хранит сами значения (ненулевые числа). Второй массив хранит индексы строк для этих значений. Третий массив хранит индексы столбцов.

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

Форматы сжатого хранения строк (CSR) и столбцов (CSC)

Индустриальным стандартом для быстрых вычислений (решения СЛАУ и умножения матрицы на вектор) стал Формат Сжатого Хранения Строк (Compressed Sparse Row, CSR). Этот формат также использует три одномерных массива, но упаковывает индексы строк гораздо умнее, чем COO. Массив значений и массив индексов столбцов остаются такими же (хранятся элементы, перебираемые слева направо и сверху вниз). Но массив строк заменяется на массив «указателей на начало строк» (row_ptr).

Массив row_ptr имеет размер, равный количеству строк в матрице плюс один. Каждый элемент этого массива показывает, на каком индексе в массиве значений начинается новая строка матрицы. Вычисляя разность между соседними элементами массива row_ptr, компьютер мгновенно узнает, сколько ненулевых элементов содержится в конкретной строке. Это обеспечивает феноменально быстрый доступ к элементам при умножении. Абсолютно симметричным является Формат Сжатого Хранения Столбцов (Compressed Sparse Column, CSC), где сжимаются индексы столбцов, а элементы записываются сверху вниз и слева направо. CSR оптимизирован для умножения y = A*x, а CSC — для y = A^T*x. Все современные библиотеки линейной алгебры (Intel MKL, SciPy.sparse, cuSPARSE на видеокартах) построены на глубокой аппаратной оптимизации этих форматов упаковки.

Подробнее

Спектральный анализ временных рядов: оконное преобразование Фурье

Проблема нестационарности в цифровой обработке сигналов

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

Однако большинство реальных сигналов в природе и технике являются нестационарными. Мелодия состоит из нот, которые звучат последовательно. Человеческая речь — это чередование согласных и гласных звуков с разным спектром. Землетрясение имеет четкое начало, фазу максимальных разрушений и затухание. Если мы применим классическое преобразование Фурье к 10-секундной аудиозаписи песни, мы получим один усредненный спектр. Этот спектр покажет нам все ноты, которые присутствовали в песне, но он абсолютно ничего не скажет о том, когда именно была сыграна каждая нота. Мы теряем информацию о времени. Для решения этой фундаментальной проблемы было изобретено Оконное преобразование Фурье (Short-Time Fourier Transform, STFT).

Скользящее окно: компромисс между временем и частотой

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

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

Весовые окна и предел Габора (принцип неопределенности)

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

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

Подробнее

Методы сглаживания данных: скользящие средние и экспоненциальное сглаживание

Борьба с высокочастотным шумом в реальных измерениях

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

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

Простое и взвешенное скользящее среднее (SMA и WMA)

Самым известным и интуитивно понятным алгоритмом является простое скользящее среднее (Simple Moving Average, SMA). Принцип его работы тривиален: мы берем «окно» заданного размера (например, 5 точек), суммируем значения внутри этого окна и делим на 5. Полученное число становится новым сглаженным значением для центральной точки окна. Затем окно сдвигается на один шаг вправо, и операция повторяется. Чем шире окно, тем сильнее сглаживается график, но тем сильнее он отстает от резких изменений реального тренда (эффект запаздывания).

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

Экспоненциальное сглаживание и фильтр Савицкого-Голея

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

Однако все методы усреднения имеют фатальный недостаток при обработке физических и химических спектров: они «срезают» острые пики полезного сигнала. Для сохранения формы пиков и точного вычисления производных в 1964 году был разработан Фильтр Савицкого-Голея. Этот метод проводит локальную полиномиальную регрессию (по методу наименьших квадратов) внутри движущегося окна. Сглаженное значение вычисляется не арифметическим усреднением, а подстановкой в найденный полином (обычно 2-й или 3-й степени). Благодаря этому фильтр Савицкого-Голея великолепно сглаживает шум, сохраняя при этом амплитуду, ширину и центры физических резонансов в спектроскопии.

Подробнее

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

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

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

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

Метод Нумерова для одномерных задач

Когда мы ищем стационарные (не зависящие от времени) состояния системы, уравнение Шредингера сводится к задаче на собственные значения. Для одномерных задач (например, частица в сложном потенциальном поле) классическим и крайне элегантным подходом является метод Нумерова. Это специализированный метод интегрирования обыкновенных дифференциальных уравнений второго порядка вида y'' = f(x)y, к которому и относится стационарное уравнение Шредингера.

Метод Нумерова базируется на разложении Тейлора и обладает впечатляющим четвертым порядком точности O(h^4), при этом требуя вычисления потенциала только в узлах сетки (без необходимости рассчитывать промежуточные точки, как в методе Рунге-Кутты). Алгоритм используется вместе с методом пристрелки (стрельбы). Мы задаем пробное значение энергии (которое является искомым собственным значением) и интегрируем уравнение с двух сторон (от плюс и минус бесконечности) навстречу друг другу. Точное значение энергии соответствует тому случаю, когда обе ветви решения гладко сшиваются в точке встречи. Этот подход до сих пор активно применяется для расчета состояний двухатомных молекул.

Временная эволюция и схема Кранка-Николсон

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

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

Подробнее

Метод характеристик для квазилинейных уравнений в частных производных

Переход от частных производных к обыкновенным

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

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

Геометрический смысл и семейство кривых

С геометрической точки зрения, решение УЧП первого порядка для функции двух переменных u(x,t) можно представить как интегральную поверхность в трехмерном пространстве (x, t, u). Метод характеристик ищет кривые, лежащие на этой поверхности. Эти кривые называются характеристиками. Вдоль характеристики изменение искомой функции происходит только за счет изменения одного параметра (часто это время t или некий фиктивный параметр s вдоль кривой).

Математически алгоритм требует составления так называемой характеристической системы ОДУ. Решая эту систему (аналитически, если это возможно, или численно, с помощью методов Рунге-Кутты), мы получаем семейство характеристических кривых на плоскости (x,t). Вдоль каждой такой кривой значение искомой функции u либо остается постоянным (для простого уравнения переноса), либо меняется по известному простому закону. Зная начальные условия задачи (профиль функции в момент времени t=0), мы можем «протащить» эти значения вдоль характеристик в любую точку будущего времени, тем самым построив полное решение исходной задачи.

Формирование ударных волн и опрокидывание

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

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

Подробнее