Задача Коши и потребность в численных методах

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

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

Метод Эйлера: простота и геометрический смысл

Исторически первым и самым простым численным методом решения ОДУ является метод Эйлера, предложенный великим Леонардом Эйлером еще в 18 веке. Этот метод основан на разложении функции в ряд Тейлора и удержании только первых двух членов разложения.

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

Семейство методов Рунге-Кутты

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

Идея методов Рунге-Кутты заключается в том, чтобы на каждом шаге интегрирования вычислить значения функции правой части (наклоны касательных) не только в начальной точке отрезка, но и в нескольких промежуточных точках (пробных точках) внутри текущего интервала. В классическом методе RK4 вычисляются четыре таких коэффициента. Итоговое приращение функции вычисляется как взвешенная сумма этих четырех коэффициентов. Метод RK4 обладает глобальной погрешностью порядка O(h^4). Это значит, что уменьшение шага всего в два раза приводит к уменьшению ошибки в 16 раз. Благодаря идеальному балансу между вычислительной сложностью и высокой точностью, метод RK4 стал стандартом де-факто («рабочей лошадкой») во многих программных пакетах инженерного анализа.

Почему прямых методов бывает недостаточно?

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

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

Метод простой итерации (Метод Якоби)

Метод Якоби является концептуально самым простым итерационным методом. Идея заключается в том, чтобы выразить из первого уравнения системы первую неизвестную x_1, из второго — x_2, и так далее до x_n. Таким образом, мы получаем систему уравнений, где слева стоят неизвестные, а справа — выражения, зависящие от всех остальных неизвестных.

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

Метод Гаусса-Зейделя и условия сходимости

Метод Гаусса-Зейделя представляет собой логическое улучшение метода Якоби. Немецкий математик Филипп Людвиг фон Зейдель заметил, что в процессе вычисления нового вектора мы можем использовать уже обновленные значения компонент. То есть, вычислив новую компоненту x_1, при вычислении x_2 мы подставляем не старое значение x_1 с предыдущей итерации, а только что найденное, более точное. Благодаря этому метод Гаусса-Зейделя сходится в среднем в два раза быстрее метода Якоби и требует меньше оперативной памяти для хранения промежуточных векторов.

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

Когда необходимо численное дифференцирование?

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

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

Формулы конечных разностей

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

• Правая (левая) разность. Вычисляется по формуле: f_prime(x) приближенно равно (f(x+h) - f(x)) / h. Это простейшая формула, имеющая порядок точности O(h). Она несимметрична и вносит определенное смещение в результат вычислений, но часто используется на границах интервалов, где данные с одной из сторон отсутствуют. Аналогично записывается левая разность: (f(x) - f(x-h)) / h.
• Центральная разность. Эта формула использует точки по обе стороны от точки интереса: f_prime(x) приближенно равно (f(x+h) - f(x-h)) / (2*h). За счет симметрии центральная разность обладает гораздо более высокой точностью — второй порядок аппроксимации O(h^2). Это означает, что при уменьшении шага h в 10 раз, ошибка метода уменьшится в 100 раз. Центральные разности являются золотым стандартом в численном моделировании при решении дифференциальных уравнений в частных производных.

Парадокс численного дифференцирования

На первый взгляд кажется, что для получения идеальной точности нужно просто взять шаг h как можно ближе к машинному нулю. Однако здесь возникает знаменитый парадокс численного дифференцирования: операция дифференцирования является некорректно поставленной задачей по Адамару.

При уменьшении шага h погрешность метода (усечения ряда Тейлора) действительно уменьшается. Но в формулах разностей в числителе мы вычитаем два очень близких числа: f(x+h) и f(x). При ограниченной разрядности ЭВМ это приводит к так называемой потере значащих цифр из-за катастрофического взаимного уничтожения. Разность близких чисел обрастает «цифровым шумом», а затем эта ошибка делится на очень малое число h, что приводит к гигантскому взрыву погрешности округления.

Таким образом, график полной погрешности имеет форму буквы «V». Слева погрешность велика из-за ошибок округления, справа — из-за математического усечения. Задача инженера или математика — найти оптимальный шаг h, находящийся на дне этой впадины, где сумма обеих ошибок минимальна. Для 64-битных систем (double precision) оптимальный шаг обычно лежит в районе квадратного корня из машинного эпсилон (примерно 10^-8).