Computer-algebraic approach to first differential approximation: Van der Pol oscillator

Аннотация

First differential approximation has been used to analyze various numerical methods for solving systems of ordinary differential equations. This has made it possible to estimate the stiffness of the ODE system that models the oscillations of the Van der Pol oscillator and the error of the method as well as to propose simple criteria for choosing a calculation step. The presented methods allow one to perform efficient calculations using computer algebra systems.

Толық мәтін

1. Введение

В 60-х годах прошлого века Н.Н. Яненко [1] сформулировал метод дифференциальных приближений разностной схемы. Основная идея этого метода состоит в замене исследования свойств разностной схемы исследованием некоторой задачи с дифференциальными уравнениями, занимающими промежуточное положение между исходной дифференциальной задачей и аппроксимирующей ее разностной схемой. В работах Н.Н. Яненко и его учеников в результате были сформулированы такие понятия, как аппроксимационная вязкость разностной схемы и первое дифференциальное приближение разностной схемы.

Первое дифференциальное приближение (ПДП) для дифференциальных уравнений в частных производных эволюционного типа, в частности, уравнения Кортевега-де Фриза с использованием систем компьютерной алгебры рассмотрено в [16], а для уравнений Навье–Стокса в [17] .

В [18] рассмотрено ПДП для разностных схем, описывающих обыкновенные дифференциальные уравнения. Обсуждена связь между сингулярным возмущением исходной системы и понятием ПДП. Для этого простого случая показана связь между методом оценки ошибки аппроксимации решения, основанным на анализе ПДП, и методом Ричардсона–Калиткина. Приведены примеры, когда при аппроксимации совместной системы дифференциальных уравнений в частных производных не всегда получаются совместные разностные системы уравнений. В качестве способа проверки совместности системы разностных уравнений предлагается проверять совместность ПДП для разностной системы. Рассмотрены вопросы вычисления ПДП в системах компьютерной алгебры Sage и SymPy.

Существует большое количество численных методов решения систем обыкновенных дифференциальных уравнений (ОДУ). Применение для их исследования ПДП позволяет получить информацию о качестве выбранного численного метода для конкретной системы, используя только символьные вычисления. В данной работе будут рассмотрены методы Рунге–Кутты [19] и многошаговые методы Адамса–Башфорта и Адамса–Мультона [20, 21] на примере осциллятора Ван дер Поля [22].

2. ПДП для систем ОДУ

Для систем ОДУ, в отличие от общего случая дифференциальных уравнений, построение ПДП не требует применения специальных алгоритмов. Обычно для численных методов требуется система ОДУ первого порядка, разрешенная относительно первых производных. Алгоритм построения ПДП для такого вида систем представляет набор простых операций с формальными степенными рядами.

Рассмотрим систему ОДУ первого порядка, разрешенную относительно производных. Также наложим условие, что правая часть уравнений F_i представляет собой полиномы от искомых функций u¹, …uⁿ:

$\left\{\begin{array}{l}{u}_t^1{F}_1\mathrm{,}\\ \dots \\ u_t^n{F}_n\mathrm{.}\end{array}\right.$ (2.1)

Применяя к системе (2.1) конкретный численный метод и заменяя искомые функции их разложением в ряд Тейлора с шагом h, получим в любой выбранной точке набор уравнений, которые будут представлять собой равенство нулю степенных рядов по h с коэффициентами G_{j, k} от искомых функций и их производных

$\left\{\begin{array}{l}{u}_t^1-F_1+\sum _{k=0}^{\infty }{G}_{\mathrm{1,}k}{h}^k=0\\ \dots \\ u_t^n-F_n+\sum _{k=0}^{\infty }{G}_{n\mathrm{,}k}{h}^k=0\end{array}\right.$ (2.2)

Система (2.2) образует дифференциальный идеал [23] относительно операций сложения, умножения и дифференцирования рядов, равных нулю. Сам вид рядов (2.2), разрешенных относительно первых производных, позволяет выполнить замену производных от функций u¹,…uⁿ в коэффициентах G_{j, k} через сами функции.

Можно сформулировать алгоритм построения ПДП, используя алгоритмы построения базисов Грёбнера [24] для рядов. Алгоритм построения ПДП оканчивается при получении первого ненулевого члена по степеням h, выраженного только через сами функции. Как следствие, вид ПДП не зависит от точки разложения в ряд Тейлора и представляет собой канонический вид.

На языке Python была написана программа с использованием системы компьютерной алгебры SymPy, которая строит ПДП для систем ОДУ, разрешенные относительно старших производных с полиномиальной правой частью. Головной модуль имеет имя fda_ode.py. Он и остальные файлы программ и файлы с приведенными ниже расчетами расположены по адресу https://github.com/blinkovua/sharing-blinkov/tree/master/FDA_ODE. Рассмотрим для примера уравнение Дуффинга [25], которое представляет собой важный пример системы (с одной степенью свободы) с нелинейной восстанавливающей силой $Fuu+\epsilon u^3$ с малым параметром e << 1.

$\left\{\begin{array}{l}{u}_{tt}+u+\epsilon u^3=0\\ u_t^2/2+u^2/2+\epsilon u^4/4=C\mathrm{.}\end{array}\right.$ (2.3)

Уравнение (2.3) имеет первый интеграл с константой C.

Для применения методов Рунге–Кутты перепишем (2.3) в виде системы первого порядка

$\left\{\begin{array}{l}{u}_t-u=0\\ u{1}_t+u+\epsilon u^3=0\\ u{1}^2/\text{\hspace{0.05em}}2+u^2/\text{\hspace{0.05em}}2+\epsilon u^4/4=\text{\hspace{0.05em}}C\mathrm{.}\end{array}\right.$ (2.4)

Пример использования приведен ниже.

В строке 1 произведен импорт основного модуля с вызовом необходимых библиотек, включая и систему компьютерной алгебры SymPy. Далее в строках 2–3 и 6 описываются основные символьные переменные и задается само уравнение Дуффинга.

В строках 8-10 производится проверка первого интеграла системы.

На строке 11 функция init первым аргументом принимает список зависимых функций от второго аргумента t с шагом h. Функция init также определяет представление мономов в виде списка степеней, в данном случае порядок производной по времени, и номера зависимой функции, что стандартно при построении базисов Грёбнера [10].

Первый аргумент всегда должен быть больше второго, который обычно должен быть больше предполагаемого порядка численного метода, а последний влияет только на объем вычислений.

В строках 13-24 задан метод Рунге–Кутты [5] для автономных систем. Правая часть для системы ОДУ задается в строках 12-17 списком аргументов y[0], y[1] = u, u1. Для значительного ускорения громоздких символьных вычислений применяется функция clip_n, которая обрезает ряд до заданной степени h.

В строках 25-31 происходит вычисление разложения метода Рунге-Кутты в ряд Тейлора в точке t = 0

$\left\{\begin{array}{l}-u_1+u_t+h\left(\frac{\epsilon u^3}{2}+\frac{u+u_{tt}}{2}\right)\dots =0\\ \epsilon u^3+u+u{1}_t+h\left(\frac{3\epsilon u^2{u}_1}{2}+\frac{u_1+u{1}_{tt}}{2}\right)\dots =0\end{array}\right.$ (2.5)

В строках 32–38 происходит приведение ПДП к каноническому виду, который, в частности, дает точный порядок численного метода:

$\left\{\begin{array}{l}{u}_t-u_1+h^4\left(\frac{{\epsilon }^2{u}^4{u}_1}{10}+\right.\\ \left.+\frac{\epsilon u_1\left(3u-u_1\right)\left(3u+u_1\right)}{120}+\frac{u_1}{120}\right)+\mathcal{O}{h}^5=0\\ u{1}_t+u+\epsilon u^3+h^4\left(-\frac{3{\epsilon }^3{u}^7}{80}-\right.\\ \left.-\dots -\frac{\epsilon u\left(5u^2+18u_1^2\right)}{240}-\frac{u}{120}\right)+\mathcal{O}{h}^5=0\end{array}\right.$ (2.6)

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

Параметр head=False показывает, что не нужно редуцировать первый член разложения. При его значении True будет происходить, если возможно, редукция головного члена. Это необходимо в том случае, если необходимо посчитать S-полином, как, например, сделано в работе [17] или посчитать первый интеграл, как показано ниже.

В строках 39–43 происходит вычисление разложения первого интеграла из системы (2.4) в ряд Тейлора в точке t = 0.

$\begin{array}{c}\epsilon u^3{u}_t+u{u}_t+u_{1}u{1}_t+h\times \\ \times \left(\begin{array}{c}\frac{\epsilon u^2\left(u{u}_{tt}+3u_t^2\right)}{2}+\\ +\frac{u{u}_{tt}+u_{1}u{1}_{tt}+u{1}_t^2+u_t^2}{2}\end{array}\right)\dots =0.\end{array}$ (2.7)

В строках 45–47 происходит приведение первого интеграла к каноническому виду. ПДП первого интеграла позволяет осуществлять дополнительный контроль численного метода для конкретного вида ОДУ (2.8).

$h^4\left(\begin{array}{c}-\frac{{\epsilon }^3{u}^7{u}_1}{16}-\frac{{\epsilon }^2{u}^3{u}_1\left(3u^2-5u_1^2\right)}{24}-\\ -\frac{\epsilon u{u}_1\left(3u^2-4u_1^2\right)}{48}\end{array}\right)+\mathcal{O}\left(h^5\right)=0$ (2.8)

3. ПДП для осциллятора Ван дер Поля

Запишем уравнение осциллятора Ван дер Поля

$u_{tt}-\mu \left(1-u^2\right)u_t+u=0$ (3.1)

в виде системы двух уравнений

$\left\{\begin{array}{l}{u}_t-u_1\\ {u_1}_t-\mu \left(1-u^2\right)u_1+u=0\end{array}\right.$ (3.2)

Применяя метод Рунге–Кутты четвертого порядка для (3.2), получим

$\left\{\begin{array}{l}{u}_1+h\left(-\frac{\mu u_1\left(u-1\right)\left(u+1\right)}{2}-\frac{u}{2}\right)+\dots =0\\ -\mu u_1\left(u-1\right)\left(u+1\right)-u+\\ +h\left(\frac{{\mu }^2{u}_1{\left(u-1\right)}^2{\left(u+1\right)}^2}{2}+\right.\\ \left.+\frac{\mu u\left(u^2-2u_1^2-1\right)}{2}-\frac{u_1}{2}\right)+\dots =0\end{array}\right.$ (3.3)

Разложение в ряд Тейлора при t = 0 для (3.3) дает

$\left\{\begin{array}{l}{u}_t-u_1+h\left(\frac{\mu u_1\left(u-1\right)\left(u+1\right)}{2}+\right.\\ \left.+\frac{u+u_{tt}}{2}\right)+\mathcal{O}{h}^2=0\\ {u_1}_t-\mu u_1\left(1-u^2\right)+u+\\ +h\left(-\frac{{\mu }^2{u}_1{\left(u-1\right)}^2{\left(u+1\right)}^2}{2}-\right.\\ \left.-\frac{\mu u\left(u^2-2u_1^2-1\right)}{2}+\frac{u_1+u{1}_{tt}}{2}\right)+\mathcal{O}\left(h^2\right)=0\end{array}\right.$ (3.4)

Используя алгоритмы построения базисов Грёбнера для рядов в системе компьютерной алгебры SymPy, построим ПДП для (3.4), в результате получим

$\left\{\begin{array}{l}{u}_t-u_1+h^4\left(\frac{{\mu }^4{u}_1{\left(u-1\right)}^4{\left(u+1\right)}^4}{120}+\right.\\ \left.+\dots \right)+\mathcal{O}\left(h^5\right)=0\\ {u_1}_t-\mu u_1\left(1-u^2\right)+u+\\ +h^4\left(-\frac{{\mu }^5{u}_1{\left(u-1\right)}^5{\left(u+1\right)}^5}{120}-\right.\\ \left.-\dots \right)+\mathcal{O}\left(h^5\right)=0\end{array}\right.$ (3.5)

Построение ПДП позволило правильно определить порядок численного метода и его остаточный член. Были построены ПДП для различных явных и неявных Рунге–Кутты и многошаговых методов Адамса–Башфорта и Адамса–Мультона. Все результаты показали, что ПДП второго уравнения системы (3.2) имеет вид $h^{p}C{\mu }^{p+1}{u}_1{u^2-}^{p+1}+\dots$, где p – порядок метода, C – некоторая константа. Значения констант C представлены в табл. 1.

Таблица 1

Метод	p	C
Гаусса–Лежандра (неявное правило средней точки)	2	1/12
Кранка–Николса [27]	2	1/12
Рунге–Кутты	4	–1/120
Адамса–Башфорта	4	–251/720
Адамса–Мультона	4	19/720
Дормана–Принса [28]	5	–1/3600
Адамса–Башфорта	5	95/288
Адамса–Мультона	5	19/720

Вычисление всех представленных примеров на компьютере с процессором Intel(R) Core(TM) i7-4770K CPU@3.50GHz в системе Debian 12 в оболочке Jupiter потребовало 620 сек и 160 Мегабайт оперативной памяти.

Остаточный член показывает жесткость системы [26] относительно параметра m. В результате необходимо выбирать шаг h таким образом, чтобы обеспечить малость остаточного члена. Вычисленный остаточный член для методов Рунге–Кутты позволяет организовать эффективное определение размера шага h для заданной погрешности.

4. Численные эксперименты

Были проведены расчеты для всех представленных примеров с начальными условиями u = 0, u_t = 1 при t = 0 до t = 100. Обозначение max показывает максимальное отклонение значения u от 0 и выбрано для контроля точности. Для табл. 2 шаг задан следующим значением h = 0.01. Символ * обозначает прерывание вычислений при значении abs(u) >= 3.

Таблица 2

μ	h pμ p+1	max	μ	h pμ p+1	max
Гаусса–Лежандра 2			Кранка–Николса 2
0	0.00e+00	1.000	0	0.00e+00	1.000
1	1.00e-04	2.009	1	1.00e-04	2.009
5	1.25e-02	2.022	5	1.25e-02	2.022
10	1.00e-01	2.015	10	1.00e-01	2.015
20	8.00e-01	2.011	20	8.00e-01	2.011
50	1.25e+01	2.018	50	1.25e+01	2.018
100	1.00e+02	*	100	1.00e+02	*
150	3.38e+02	*	150	3.38e+02	*
Рунге–Кутты 4			Дормана–Принса 5
0	0.00e+00	1.000	0	0.00e+00	1.000
1	1.00e-08	2.009	1	1.00e-10	2.009
5	3.13e-05	2.022	5	1.56e-06	2.022
10	1.00e-03	2.014	10	1.00e-04	2.014
20	3.20e-02	2.008	20	6.40e-03	2.008
50	3.12e+00	2.004	50	1.56e+00	2.001
100	1.00e+02	*	100	1.00e+02	2.265
150	7.59e+02	*	150	1.14e+03	*
Адамса–Башфорта 4			Адамса–Мультона 4
0	0.00e+00	1.000	0	0.00e+00	1.000
1	1.00e-08	2.009	1	1.00e-08	2.009
5	3.13e-05	2.022	5	3.13e-05	2.022
10	1.00e-03	2.014	10	1.00e-03	2.014
20	3.20e-02	*	20	3.20e-02	2.008
50	3.12e+00	*	50	3.12e+00	2.008
100	1.00e+02	*	100	1.00e+02	*
150	7.59e+02	*	150	7.59e+02	*
Адамса–Башфорта 5			Адамса–Мультона 5
0	0.00e+00	1.000	0	0.00e+00	1.000
1	1.00e-10	2.009	1	1.00e-10	2.009
5	1.56e-06	2.022	5	1.56e-06	2.022
10	1.00e-04	*	10	1.00e-04	2.014
20	6.40e-03	*	20	6.40e-03	2.008
50	1.56e+00	*	50	1.56e+00	2.022
100	1.00e+02	*	100	1.00e+02	*
150	1.14e+03	*	150	1.14e+03	*

Как видно из табл. 2, значения h ^pμ^p+1 может использоваться для выбора h при вычислениях с фиксированным шагом. Также заметно преимущество неявных методов, по сравнению с явными методами, что возможно обусловлено меньшим значением константы C в ПДП.

Для вычислений с переменным шагом для заданной погрешности tol = 0.01 по следующей формуле $\min htol+abs{\text{остаточный член}}^p$. Как показано в табл. 3, при начальном шаге h = 0.01 все вычисления показали устойчивость при изменении шага и большую сходимость при большем порядке метода.

Таблица 3

μ	шагов	max	μ	шагов	max
Гаусса–Лежандра 2			Кранка–Николса 2
0	10001	1.000	0	10001	1.000
1	10001	2.009	1	10001	2.009
5	12999	2.022	5	14079	2.022
10	14922	2.014	10	16211	2.014
20	16783	2.008	20	18441	2.008
50	16899	2.003	50	18471	2.003
100	20332	2.001	100	22327	2.001
150	32912	2.001	150	35752	2.001
Рунге–Кутты 4			Дормана–Принса 5
0	10001	1.000	0	10001	1.000
1	10001	2.009	1	10001	2.009
5	10001	2.022	5	10001	2.022
10	10024	2.014	10	10001	2.014
20	10179	2.008	20	10004	2.008
50	10225	2.003	50	10041	2.003
100	10273	2.001	100	10062	2.001
150	11796	2.001	150	10815	2.001

5. Заключение

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

Применяя ПДП, возможно обнаружить и оценить жесткость системы ОДУ, а также использовать остаточный член ПДП для выбора постоянного шага или для управления переменным шагом.

Исходные тексты программ и представленные вычисления приведены по адресу github.com/blinkovua/sharing-blinkov/tree/master/FDA_ODE .

Представленные методы позволяют провести эффективные вычисления средствами систем компьютерной алгебры для решений систем ОДУ с правой полиномиальной частью.