Смешанная символьно-числовая теория возмущений

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

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

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

Отличительной чертой метода возмущений является генерация длинных и запутанных промежуточных уравнений, которые подвергаются алгоритмическим и механическим операциям. Поэтому такие задачи хорошо подходят для CAS. Фактически программные пакеты CAS используются для помощи в расчетах возмущений с начала 1970-х годов.

В этом руководстве мы хотим показать, как использовать сочетание символьных операций (Symbolics.jl) и численных методов для решения простых задач возмущения.

Начнем с аналога "hello world!" в задачах теории возмущений, решая уравнения пятого порядка. Нам нужно найти вещественное значение , такое, что . Согласно теореме Абеля, общее уравнение пятого порядка не имеет замкнутого решения. Конечно, мы можем легко решить это уравнение численно, например с помощью метода Ньютона. Используем следующую реализацию метода Ньютона.

В этом коде Symbolics.derivative(eq, x) делает именно то, что следует из его названия: вычисляет символьную производную eq (выражение Symbolics.jl) по x (переменной Symbolics.jl). Мы используем Symbolics.substitute(eq, D) для определения формулы обновления, подставляя переменные или подвыражения (определенные в словаре D) в eq. Следует отметить, что substitute является исполнительным компонентом нашего кода и будет использоваться много раз в остальных частях этого руководства. solve_newton написан с учетом простоты и ясности, а не производительности.

Вернемся к уравнению пятого порядка. Мы можем определить переменную Symbolics как @variables x, а затем решить уравнение solve_newton(x^5 + x - 1, x, 1.0) (здесь x₀ = 0 — наше первое предположение). Ответ: 0,7549. Теперь посмотрим, как можно решить ту же задачу с помощью методов возмущения.

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

является решением простого уравнения, поэтому . Подставляем в исходную задачу,

Разложив уравнения, получаем

Это уравнение должно выполняться для каждой степени . Таким образом,

и

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

И, наконец,

При решении исходной задачи получаем , по сравнению с 0,7548, рассчитанным численно. Мы можем повысить точность, включив больше членов в разложение . Однако вычисления, хотя и являются простыми, очень быстро становятся запутанными и трудновыполнимыми для проведения вручную. Вот почему система CAS очень полезна для решения задач возмущения.

Теперь посмотрим, как можно выполнить те же вычисления в Julia. Пусть будет порядком разложения. Начнем с определения символьных переменных:

Затем определим

Следующим шагом будет подстановка x в уравнение задачи

Разложенная форма eq имеет вид

Нам нужен способ получения показателей различных степеней ϵ. Функция collect_powers(eq, x, ns) возвращает степени переменной x в выражении eq. Аргумент ns представляет собой диапазон степеней.

Чтобы вернуть показатели и в eq, мы можем написать

Несколько слов о том, как работает функция collect_powers. Она использует substitute для нахождения показателя заданной степени x, передавая Dict со всеми степенями x, которым задано значение 0, кроме целевой степени, которой задано значение 1. Например, следующее выражение возвращает показатель ϵ^2 в eq,

Вернемся к задаче. При наличии показателей степеней ϵ мы можем задать каждое уравнение в eqs равным 0 (помните, мы перестроили задачу так, что eq равен 0) и решить систему линейных уравнений, чтобы найти числовые значения показателей. В Symbolics.jl есть функция Symbolics.solve_for, которая может решать системы линейных уравнений. Однако присутствие членов более высокого порядка в eqs препятствует корректной работе Symbolics.solve_for(eqs .~ 0, a). Мы можем воспользоваться тем фактом, что система имеет треугольную форму, и начать с решения eqs[1] при a₁, затем подставить это в eqs[2] и решить при a₂ и так далее. Этот каскадный процесс реализуется с помощью функции solve_coef(eqs, ps):

Здесь eqs является массивом выражений (предполагается, что они равны 0), а ps — массивом переменных. В результате получаем словарь пар переменная => значение. Применим solve_coef к eqs, чтобы получить численные значения параметров:

Наконец, подставим обратно значения a в определение x как функцию от 𝜀. Обратите внимание, что 𝜀 — это число (обычно Float64), а ϵ — символьная переменная.

Таким образом, решением нашей исходной задачи становится X(1), что равно 0,76. Можно использовать большие значения n, чтобы повысить точность оценок.

Помните, что числовое значение равно 0,7549. Во всех примерах этого и следующего руководства используются две функции — collect_powers и solve_coef(eqs, a).

Исторически методы возмущений были впервые изобретены для выполнения орбитальных расчетов для Луны и планет. Отдавая дань этой истории, наш второй пример связан с космосом. Наша цель — решить уравнение Кеплера:

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

Как и в первом примере, эту задачу легко решить с помощью метода Ньютона. Например, пусть (эксцентриситет Земли) и . solve_newton(x - e*sin(x) - M, x, M) равен 1,5875 (по сравнению с π/2 = 1,5708). Теперь попробуем решить ту же задачу с помощью методов возмущения (см. функцию test_kepler).

При получаем . Поэтому можем использовать в качестве параметра возмущения. Для согласованности с другими задачами также переименуем в , а в .

С этого момента мы используем вспомогательную функцию def_taylor для определения ряда Тейлора, вызывая ее как x = def_taylor(ϵ, a, 1), где аргументами являются, соответственно, переменная возмущения, которая представляет собой массив показателей (начиная с показателя ), и необязательный постоянный член.

Начнем с определения символьных переменных (предполагая, что n = 3):

Еще больше упрощаем задачу, заменяя sin степенным рядом с помощью вспомогательной функции expand_sin:

Чтобы протестировать,

Уравнение задачи имеет следующий вид:

Действуем так же, как и в первом примере. Собираем показатели степеней ϵ

и затем решаем с a:

Наконец, подставляем vals обратно в x:

Результат — 1,5876, по сравнению с числовым значением 1,5875. Принято задавать порядок X на основе степеней 𝑀, а не 𝜀. Мы можем вычислить этот ряд как collect_powers(sol, M, 0:3). Результат (после очистки вручную) таков:

Сравнение формулы с формулой для xD835__xDC38 в статье Википедии об уравнении Кеплера:

Первое отклонение — в показателе .