CC BY
4
УДК 519.62
БОТ: 10.18698/2309-3684-2022-3-1832
Применение одношагового метода Галеркина для решения системы обыкновенных дифференциальных уравнений с начальными условиями
© СВ. Русских, Ф.Н. Шклярчук Московский авиационный институт, Москва, 125993, Россия
Рассматривается нелинейная колебательная система, описываемая обыкновенными дифференциальными уравнениями с переменными коэффициентами. Предполагается, что на рассматриваемом интервале времени решение системы является достаточно гладкими — без разрывов, столкновений и бифуркаций. Из неоднородной системы уравнений выделяются в явном виде члены, линейно зависящие от координат, скоростей и ускорений и члены, зависящие от этих переменных нелинейно. Предлагается новый подход для численного решения шаговым методом начальной задачи, описываемой такой системой обыкновенных дифференциальных уравнений второго порядка. На шаге интегрирования неизвестные функции представляются в виде суммы функций, удовлетворяющих начальным условиям: линейного решения Эйлера и нескольких заданных корректирующих функций в виде полиномов второй и выше степеней с неизвестными коэффициентами. Дифференциальные уравнения на шаге удовлетворяются приближённо в смысле слабого решения по методу Галеркина на системе корректирующих функций. Получаются алгебраические уравнения с нелинейными членами, которые решаются методом итераций, начиная в первом приближении с линейного решения. Полученное решение в конце данного шага используется в качестве начальных условий на последующем шаге. В качестве примера рассмотрено одно однородное дифференциальное уравнение второго порядка без первой производной с сильной кубической нелинейностью по координате (при максимальной амплитуде нелинейная сила в два раза превышает линейную силу). Это уравнение имеет точное периодическое решение в виде интеграла энергии консервативной системы, которое используется для оценки точности численных решений, полученных методами Галеркина, Рунге-Кутта и Адамса второго порядка, а также методами Radau5 и BDF на различных интервалах времени (до 8000 периодов свободных колебаний системы) при использовании различных постоянных шагов интегрирования (от 0,0025 долей периода). При этом в методе Галеркина на каждом шаге использовалось четыре одинаковых корректирующих функций в виде полиномов от второй до пятой степеней. Показано, что на больших интервалах времени вычислений метод Галеркина обладает более высокой точностью по сравнению с другими рассмотренными численными методами. Поэтому он может быть использован для численного решения нелинейных задач, в которых требуется решать их на больших интервалах времени; например при расчете установившихся предельных циклов нелинейных колебаний и хаотических нелинейных колебаний со странными аттракторами.
Ключевые слова: обыкновенные дифференциальные уравнения, нелинейные системы, начальная задача, численные решения, одношаговый метод Галеркина
Введение. Данная Численным методам решения обыкновенных дифференциальных уравнений (ОДУ) посвящена обширная литература: учебные пособия [1], [2] и др.; монографии [3-5] и др.;
многочисленные научные статьи. В фундаментальных трудах [4], [5] на основе исследований авторов и других исследователей (приведены ссылки на сотни публикаций) представлены теоретические обоснования вычислительных схем различных модификаций и комбинаций явных и неявных одношаговых методов типа Рунге-Кутта и многошаговых методов типа Адамса и некоторых других методов для нежестких и жестких задач, включая дифференциально-алгебраические задачи. Большое внимание уделено выбору шага, оценкам локальных и глобальных (накопленных) погрешностей, а также оценкам устойчивости и сходимости численных решений. Приведено много примеров решения тестовых задач с использованием различных компьютерных программ со сравнениями и анализом результатов. В этих работах в основном рассматриваются системы первого и второго порядков, разрешенные относительно старших производных. В одном примере [5] для многозвенной механической системы с матрицей инерции, зависящей от времени и перемещений системы, ускорения рассматриваются как независимые параметры и уравнения движения системы записываются в виде нелинейных алгебраических уравнений, к которой добавляются дифференциальные соотношения для скоростей и перемещений. Наряду с многочисленными работами, отмеченными в [5], жестким и дифференциально-алгебраическим задачам также посвящены работы [6-10] и др.
Отметим некоторые оригинальные подходы к численному решению начальной задачи для ОДУ. В работе [11] предлагается парамет-ризировать уравнения и затем интегрировать их по параметру. Метод наименьших квадратов для минимизации невязки на шаге применен к решению линейных задач для ОДУ и для дифференциально-алгебраических задач с использованием разложения в ряд по заданным функциям с неизвестными коэффициентами в работах [12], [13]. В работах [14], [15] и ряде других работ этих авторов для приближенного решения начальной задачи, описываемой ОДУ, используются разложения по полиномам Чебышева, а в работах [16], [17] — по полиномам Эрмита.
Метод Галеркина обычно используется при решении краевых задач в слабой формулировке для сведения их к алгебраическим уравнениям или к ОДУ на системе глобальных или локальных базисных функций [1-3]; в последнем случае этот метод также как в вариационной формулировке задачи называется методом конечных элементов или методом конечных объемов [18], [19].
Большие нестационарные системы, в которых одновременно происходят медленные и быстрые движения, с одной стороны требуют выполнения расчетов на больших интервалах времени, а с другой стороны — использования весьма малых шагов интегрирования по времени (меньше периода колебаний системы с максимальной собственной частотой). Например, в случае дискретных моделей
неоднородных сплошных сред с отраженными ударными волнами, локальными разрушениями и пр. для получения достоверного решения необходимо выполнить численные эксперименты, используя различные методы (программы) и вычисления на разных шагах. Это также относится к нелинейным системам при расчете предельных циклов колебаний, а также хаотических колебаний со странными аттракторами
Вопросы сходимости численных решений ОДУ и их локальной и глобальной погрешностей, а также выбора длины шага исследовались в работах [4], [5], отмеченных там в соответствующих ссылках работах и в [10], [20], [21]. При решении ОДУ большой размерности, записанных в общей матричной форме, которые, например, получаются после дискретизации по координатам в задачах динамики сплошных деформируемых сред и тел (конструкций) и нестационарной аэрогидродинамики [18], [22], для дискретизации по времени обычно используются сравнительно простые схемы: 1) центральные разности второго порядка с квадратичной аппроксимацией на двух шагах [23]; 2) нецентральные разности назад третьего порядка с кубической аппроксимацией на трех шагах [24]; 3) метод обобщенного ускорения с квадратичной аппроксимацией на одном шаге и с двумя выбираемыми из условия устойчивости параметрами [25]; 4) модифицированный метод линейного ускорения (9 -метод), в котором считается, что ускорение изменяется линейно на расширенном шаге 9Н, где 9> 1,37 для безусловной устойчивости метода [26]. Эти четыре метода приводят задачу к рекуррентным алгебраическим соотношениям в матричном виде, которые решаются с учетом начальных условий по шагам; они являются частными случаями общей разностной схемы [27].
Постановка задачи и метод её решения. Рассмотрим нелинейную систему обыкновенных дифференциальных уравнений с переменными коэффициентами, приведенную к виду
Ь(х,х,х,/) = тх + с!х + сх + ]Ч-Х = 0, (1)
где х(/) = {хг} — вектор неизвестных; Х(/) = {Хг} — вектор правых частей неоднородной системы дифференциальных уравнений; ^х, X, X, t) = {Ni } — вектор нелинейных членов произвольного вида, который определяется конкретной рассматриваемой задачей; т = [ту (0], й = [йу (0], с = [е^ (/)] — квадратные матрицы инерции,
демпфирования и жесткости соответственно. Здесь и далее точкой обозначается производная по времени t.
Начальные условия для вектора неизвестных координат и их скоростей в момент времени t = ^ задаются в виде:
х(0 = хо, уМ = ±М = ±о- (2)
Время г рассматривается в виде возрастающей последовательности дискретных моментов гк, к — 0,1,2,___, п в общем случае с переменным шагом к — гк — гкЛ. Для удобства вычислений на каждом шаге вводим безразмерное время %:
% — 1 (г—гк_1), 0 <%< 1.
Далее штрихом обозначается производная по безразмерному времени %. Имеет место соотношение: ё — 1.
ёг к
Векторы перемещений х и скоростей V — Х на к -ом шаге аппроксимируются как
£
х — Хк—1 + Vk—-к% + Х Ж (гк ф (%);
(3)
к— 1 г ' г / К—1 1 ¿^^
Г=1 " Г=1
где первые два слагаемых в выражении для вектора х (соответственно, первое слагаемое для V) представляют явное решение Эйлера, удовлетворяющее начальным условиям (2), а заданные корректирующие функции (г (%), г — 1,2, _, £, удовлетворяющие начальным условиям (г (0) — 0, (р'г (0) — 0, берутся одинаковыми на всех шагах к — 1,2, _, п и для всех переменных х; ж(к) — неизвестные векторы коэффициентов.
Уравнения (1) на к -ом шаге (гк-1 < г < гк ) будем удовлетворять приближенно в смысле слабого решения по методу Галеркина:
1
= ¿ = 1,2,..„я, р = 1,2,...,5. (4)
о
В результате получаем систему алгебраических уравнений на к -ом шаге для векторов неизвестных коэффициентов ж(к), к —1,2, _,п,
Г — 1,2, _, £ .
В случае нелинейных уравнений система (1) в виде (4) решается по методу последовательных приближений на каждом шаге:
х^х(й), >7 = 1,2,...:
в качестве первого приближения используется решение линейной задачи (п = 1, х ^ х(1), N ^ = 0 ) или линеаризованные по Е значения вектора нелинейных частей. Во втором приближении (п = 2, х —> х(2)) нелинейные члены определяются по результатам первого приближения Г^1' —» |. Процесс повторяется
до тех пор, пока не будет достигнута сходимость с заданной точностью. Такой подход особенно удобен для физически (не численно) устойчивых слабонелинейных систем с гладкими решениями, для которых процесс итераций сходится достаточно быстро.
Уравнения (4) с учетом (1) и (3) на каждом шаге к = 1,2,...,п можно записать в развернутом виде:
Я 1
X«й^+о^+ЧЧ-, Г,р = 1,2,...,, (5)
г=1
где а Гкр = [«г., р ], ар} = [ст£р], т (гк} = [т^] — определяемые на каждом шаге интегрирования матрицы коэффициентов:
а« = \(шфг+Афг+сфг)фр^ =
(6)
ар' = } еда, тГк' = } (d + ек£)Ф&.
Если функции т (Е), ^ (Е), сг>. (Е) являются достаточно
гладкими и слабо изменяются в пределах шага к, то их на каждом шаге можно считать постоянными.
После решения уравнений (5) с заданной точностью вычисления
коэффициентов г(к> затем с использованием выражений (3)
определяются вектора х и V при Е = 1 для момента времени £ = ^ :
Хк = Хк_1 + Vk+ £ X ^фг (1);
(7)
V* = V,, +£ж{?ФЛ 1) = + {±$4 (1).
г=1 " г=1
Далее полученные значения векторов хк и vk используются в качестве начальных условий для решения уравнений (5) на (к +1) -ом шаге.
Выбор корректирующих функций. В качестве корректирующих функций ф (Е), уточняющих линейный закон изменения
о
переменных на малом временном шаге, в рядах (3) можно использовать следующие наборы пяти ( г, р —1,2, _ £ — 5 ) линейно-независимых степенных функций, удовлетворяющих условиям фг (0) — 0, фф (0) — 0 .
1) Базовые степенные функции:
ф —%2, ф —%3, фз —%4, ф4 —%5, ф5 —%6 . (8)
2) Ортогональные полиномы, образованные из базовых степенных функций (8):
ф — %2;
ф2 — 6%3 — 5%2; ф— 28%4 — 42%3 +15%2; (9)
ф —120%5 — 252%4 +168%3 — 35%2; ф — 495%б —1320%5 +1260%4 — 504%3 + 70%2.
Они удовлетворяют условиям ( (1) — 1, а также условиям
1
\ффМ = 0 при р
ф г,
1 ^^
2 г + 3
3) Степенные полиномы, удовлетворяющие определенным условиям при % — 1, упрощающие спряжения решений на соседних шагах (7):
ф — 3%2 — 2%3, ф(1) — 1, ф/(1) — 0; ф2 —%3 —%2, ф2(1) — 0, ф (1) — 1; фз —%2 — 2%3 +%4, фз(1) — ф3(1) — 0; (10)
ф4 —%2 — 4%3 + 5%4 — 2%5, ф4(1) — ф4(1) — 0; ф5 — %2 — 8%3 +19%4 —18%5 + 6%б, ф5(1) — ф5(1) — 0.
Функции (9) и (10) линейно выражаются через базовые функции (8) с соответствующими матрицами преобразования. Поэтому при использовании функций (9) и (10) интегралы (6) для упрощения следует вычислить сначала для функций (8) и затем преобразовать их с соответствующими матрицами:
Ф° — Тф1), ф1) — [%2 %3 %4 %5 %6], фг) — [ф ф2) ф3) ф4 ф5)];
0
T при i = 2,3:
1 0 0 0 0 " " 3 -2 0 0 0
-5 6 0 0 0 -1 1 0 0 0
15 -42 28 0 0 T = 1 -2 1 0 0
-35 168 -252 120 0 1 -4 5 -2 0
70 -504 1260 -1320 495 1 -8 19 -18 6
Базовые функции (8) при большом их числе (£ > 8 ) приводят к плохобусловленным матрицам коэффициентов а®, поэтому эти
функции можно использовать только при малых £ . Ортогональные полиномы (9) имеют большие коэффициенты и являются «колебательными» функциями на интервале 0 < % < 1. Поэтому при вычислениях интегралов (6) с переменными коэффициентами при большом числе £ они могут привести к численным погрешностям, например, за счет малых разностей больших чисел.
Вычислительные эксперименты показывают, что при численном интегрировании на малых шагах к сходимость с достаточно высокой точностью и с весьма близкими результатами обеспечивается при использовании четырех ( £ — 4 ) корректирующих функций (9)—(10).
Примеры расчета со сравнениями. В качестве примера рассмотрим нелинейную консервативную систему с одной степенью свободы ( п — 1), свободные колебания которой описываются дифференциальным уравнением с начальными условиями
х + 100х + 200х3 =0, х(0) = 1, х(0) = 0.
Полная энергия этой системы: Е —1/2X2 + 50(х2 + х4) —100 . Отсюда получается фазовая характеристика и период колебаний:
х = ±10-^2-х2 -х4
4 г ёх 4 ёу
ш1 Д 3 "Г ш 1
4 .
10 W2 - х2 - х4 10 J0 ^2 + sin
= 0,40043095.
sin у
условие периодичности
х
(t + mT ) = х (t) при
Выполняется m = 1,2,3, ....
Для сравнения приведем результаты численного решения данной задачи с постоянным шагом h на больших интервалах времени для моментов t = 100 и t = 1000 по стандартным программам в PTC MathCad, реализующим метод Адамса (А) четвертого порядка и метод Рунге-Кутта (Р.-К.) четвертого порядка, а также — по предла-
гаемому в данной статье методу Галеркина с использованием 5 = 4 корректирующих функций в виде трех различных наборов: 1) — (8), Г1; 2) — (9), Г2; 3) — (10), Гэ.
В таблице 1 приведены результаты вычисления Е, х и X с шагом И = 0,005 , а в табл. 2 — с шагом И = 0,001. При этом во всех случаях для решения нелинейных уравнений (5) с коэффициентами (6) с точностью до 10-6 требовалось не более четырех итераций. Как видно, результаты вычислений х и X по методу Адамса при £ = 1000 сильно расходятся с результатами вычислений по методу Рунге-Кутта и по методу Галеркина; последние два метода при И = 0,001, £ < 8000 дают близкие результаты. При увеличении интервала вычислений до £ > 8000 метод Рунге-Кутта при И = 0,001 также начинает давать ошибочные результаты.
Таблица 1
Результаты вычислений при шаге к = 0,005
Метод £ = 100 £ = 1000
Е X X Е X X
А 100,0285 -0,0686 14,1274 100,2706 0,9819 3,3402
Р.-К. 99,9746 -0,1441 14,0652 99,7474 -0,3454 13,6435
Г1 100,0016 -0,1057 14,1023 100,0158 -0,5387 -12,7509
Г2 100,0004 -0,1073 14,1010 100,0042 -0,3889 -13,5128
Гэ 100,0003 -0,1074 14,1009 100,0030 -0,3735 -13,5687
Radau5 - - - - - -
ББР 100,0575 -0,0278 14,1435 100,4619 -0,7599 10,4807
Таблица 2 Результаты вычислений при шаге к = 0,001
Метод £ = 100 £ = 1000 £ = 8000
Е X Е X Е X
А 100,0278 -0,0699 100,2747 0,9839 101,1724 -0,1051
Р.-К. 100,0000 -0,1078 99,9999 -0,3322 99,9993 -0,9672
Г1 99,9999 -0,1079 99,9994 -0,3252 - -
Г2 100,0000 -0,1078 99,9998 -0,3301 - -
Гэ 100,0000 -0,1078 100,0000 -0,3333 100,0000 -0,9859
Radau5 99,7810 -0,4164 98,2622 0,7068 - -
ББР 100,0572 -0,0286 100,4680 -0,6917 - -
На рис. 1 и рис. 2 приведены результаты вычислений х(£) с шагом И = 0,001 в пределах примерно одного периода колебаний вблизи моментов времени £ = 16000 и £ = 32000 по методам Адамса (кривые 1), Рунге-Кутта (кривые 2) и Галеркина (кривые 3). Результаты метода Галеркина в пределах масштаба графиков совпадают с численным решением задачи, полученным на первом периоде при
очень мелком (h = 10 7) шаге, которое можно считать практически точным.
х (t) 1,2 0,6 0
-0,6 -1,2
16000,0 16000,1 16000,2 16000,3 t Рис. 1. Результаты вычислений с шагом h = 0,001 вблизи времени t = 16000
х (t) 1,2
0,6 0
-0,6 -1,2
32000,0 32000Д 32000,2 32000,3 / Рис. 2. Результаты вычислений с шагом h = 0,001 вблизи времени t = 32000
Накопление глобальной ошибки вычислений по методам Адамса и Рунге-Кутта происходит в основном за счет сдвига по фазе кривых, представляющих колебания с медленным изменением периода T . При этом полная энергия E и амплитуды х и х изменяются незначительно.
Численные решения рассмотренного выше дифференциального уравнения с кубической нелинейностью также были получены с помощью стандартных программ в PTC MathCad по методам Radau5 (модифицированный метод Рунге-Кутта) и BDF (метод обратного дифференцирования), предназначенных для жестких и дифференциально-алгебраических уравнений. Результаты расчетов приведены в таблице 1 и таблице 2 в двух нижних строчках. Эти методы обеспечили высокую точность только на сравнительно малых интер-
валах времени (< 50 периодов колебаний при решении с шагом h = 0,005 и < 100 периодов колебаний при шаге h = 0,001). Следует отменить, что при шаге h = 0,005 решение по методу Radau5 не сходилось из-за нарушения вычислительной точности. На больших интервалах времени решения расходятся даже на шаге h = 0,001.
Заключение. В работе рассмотрены системно-динамическое и агентное описание модели конкуренции, осуществлена попытка переноса прогнозов поведения непрерывной модели на дискретную модель. Проведённые имитационные эксперименты как показывают возможность применения прогнозов, полученных в одной модели, к другой, так и демонстрируют наличие ограничений при подобных переходах.
Предложен новый подход для численного решения начальной задачи, описываемой нелинейной системой обыкновенных дифференциальных уравнений второго порядка с переменными коэффициентами. На шаге интегрирования каждая искомая функция представляется в виде суммы явного линейного решения Эйлера, удовлетворяющего начальным условиям, и конечного ряда корректирующих степенных функций или образованных из них полиномов с неизвестными коэффициентами. Методом Галеркина получаются алгебраические уравнения для определения неизвестных коэффициентов на данном шаге с использованием в случае нелинейной системы процедуры итераций.
В качестве примера рассмотрено однородное дифференциальное уравнение с кубической нелинейностью по координате, которое при заданных начальных условиях имеет периодическое решение. Получены численные решения этой консервативной задачи на больших интервалах времени (до 80000 периодов колебаний системы) с использованием различных наборов корректирующих полиномов от второй до пятой степени включительно и двух различных постоянных шагов, приблизительно равных 0,0125 и 0,0025 долей периода.
Выполнены сравнения с численными решениями, полученными по стандартным программам в PTC MathCad, методами Адамса и Рунге-Кутта четвертого порядка на тех же самых постоянных шагах, а также методами Radau5 и BDF.
Показано, что методом Галеркина задача решается на больших интервалах (включая 80000 периодов с шагом 0.0025 периода) практически точно. Методы Адамса и Рунге-Кутта в данном случае на интервалах более 8000 периодов дают неправильные результаты в основном за счет медленного уменьшения периода колебаний и сдвига по фазе при сохранении достаточно высокой точности амплитуд. Методы Radau5 и BDF обеспечивают высокую точность только на сравнительно малых интервалах времени, а на больших интервалах их численные решения расходятся.
Работа выполнена при финансовой поддержке РНФ (проект № 22-2901206).
ЛИТЕРАТУРА
[1] Демидович Б.П., Марон И.А., Шувалова Э.З. Численные методы анализа. Москва, Наука, 1967, 368 с.
[2] Бахвалов И.В., Жидков Н.П., Кобельков Г.М. Численные методы. Москва, Лаборатория базовых знаний, 2000, 630 с.
[3] Коллати Л. Численные методы решения дифференциальных уравнений. Москва, Издательство иностранной литературы, 1953, 460 с.
[4] Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. Москва, Мир, 1990, 512 с.
[5] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. Москва, Мир, 1999, 685 с.
[6] Скворцов Л.М. Явный многошаговый метод численного решения жестких дифференциальных уравнений. Журнал вычислительной математики и математической физики, 2007, т. 47, № 6, с. 259-267.
[7] Латыпов А.Ф., Никуличев Ю.В. Численные методы решения задач Коши для жестких систем обыкновенных дифференциальных уравнений на основе многозвенных интерполяционных полиномов Эрмита. Журнал вычислительной математики и математической физики, 2007, т. 47, № 2, с. 234-244.
[8] Булатов М.В., Тыглиян А.В., Филлипов С.С. Об одном классе одношаго-вых одностадийных методов для жестких систем обыкновенных дифференциальных уравнений. Журнал вычислительной математики и математической физики, 2011, т. 51, № 7, с. 1251-1265.
[9] Скворцов Л.М. Неявный метод пятого порядка для численного решения дифференциально-алгебраических уравнений. Журнал вычислительной математики и математической физики, 2015, т. 55, № 6, с. 978-984.
[10] Белов А.А., Калиткин Н.Н. Выбор шага по кривизне для жестких задач Коши. Математическое моделирование, 2016, т. 28, № 11, с. 97-112.
[11] Волков-Богородский Д.Б., Данилин А.Н., Кузнецов Е.Б., Шалашилин В.Н. О неявных методах интегрирования начальных задач для параметризиро-ванных систем обыкновенных дифференциальных уравнений второго порядка. Журнал вычислительной математики и математической физики, 2003, т. 43, № 11, с. 1684-1696.
[12] Булатов М.В., Горбунов В.К., Мартыненко Ю.В., Конг Н.Д. Вариационные методы к численному решению дифференциально-алгебраических уравнений. Вычислительная техника, 2010, т. 15, № 5, с. 3-13.
[13] Чистяков В.Ф., Чистякова Е.В. Применение метода наименьших квадратов для решения линейных дифференциально-алгебраических уравнений. Сибирский журнал вычислительной математики, 2013, т. 16, № 1, с. 81-95.
[14] Залеткин С.Ф. Численное интегрирование обыкновенных дифференциальных уравнений с использованием ортогональных разложений. Математическое моделирование, 2010, т. 22, № 1, с. 69-85.
[15] Арушанян О.Б., Залеткин С.Ф. Об одном аналитическом методе приближенного решения канонических систем обыкновенных дифференциальных уравнений второго порядка. Вестник Московского Университета. Серия 1. Математика. Механика, 2019, № 3, с. 65-69.
[16] Аульченко С.М., Латыпов А.Ф., Никульчев Ю.В. Метод численного интегрирования систем обыкновенных дифференциальных уравнений с
использованием интерполяционных полиномов Эрмита. Журнал вычислительной математики и математической физики, 1998, т. 38, № 10, с. 1665-1670.
[17] Латыпов А.Ф., Попик О.В. Численный метод решения задачи Коши для жестких обыкновенных дифференциальных уравнений на основе многозвенных интерполяционных полиномов Эрмита. Вычислительная техника, 2011, т. 16, № 2, с.78-85.
[18] Ершов Н.Ф., Шахверди Г.Г. Метод конечных элементов в задачах гидродинамики и гидроупругости. Ленинград, Судостроение. 1984, 240 с.
[19] Ioriatti M., Dumbser M. Semi-implicit staggered discontinuous Galerkin schemes for axially symmetric viscous compressible flows in elastic tube. Computers and Fluids, 2018, vol. 167, pp. 166-179.
[20] Куликов Г.Ю., Хрусталёва Е.Ю. Об автоматическом управлении длиной шага и порядком в одношаговых коллокационных методах со старшими производными. Журнал вычислительной математики и математической физики, 2010, т. 50, № 6, с. 1060-1077.
[21] Вайнер Р., Куликов Г.Ю. Эффективное управление точностью численного интегрирования обыкновенных дифференциальных уравнений и оптимальные интерполяционные равнозначные блочные методы с переменным шагом. Журнал вычислительной математики и математической физики, 2014, т. 54, № 4, с. 591-607.
[22] Баженов В.Г., Чекмарев Д.Т. Решение задач нестационарной динамики пластин и оболочек вариационно-разностным методом. Нижний Новгород, Издательство ННГУ, 2000, 118 с.
[23] Belytschko T. A survey of numerical methods and computer programs for dynamic structural analysis. Nuclear Engineering and Design, 1976, vol. 37, iss. 1, pp. 23-34.
[24] Houbolt J.C. A recurrence matrix solution for the dynamic response of elastic aircraft. Journal of Aeronautical Sciences, 1950, vol. 17, pp. 540-550.
[25] Newmark N.M. A method of computation for structural dynamics. ASCE Journal of the Engineering Mechanics Division, 1959, vol. 85, pp. 67-94.
[26] Bathe K.J., Wilson E.L. Numerical methods in finite element analysis. New York, Prentice-Hall, 1976, 544 p.
[27] Krieg R.D., Key S.W. Transient shell response by numerical time integration. International Journal for Numerical Methods in Engineering, 1973, vol. 7, no. 3, pp. 273-286.
Статья поступила в редакцию 13.08.2022
Ссылку на эту статью просим оформлять следующим образом: Русских С.В., Шклярчук Ф.Н. Применение одношагового метода Галеркина для решения системы обыкновенных дифференциальных уравнений с начальными условиями. Математическое моделирование и численные методы, 2022, № 3, с. 18-32.
Русских Сергей Владимирович — д-р физ.-мат. наук, доцент, доцент кафедры «Проектирование и прочность авиационно-ракетных и космических изделий», Московского авиационного института. e-mail: sergey.russkih@rambler.ru
Шклярчук Федор Николаевич — д-р. техн. наук, профессор кафедры «Проектирование и прочность авиационно-ракетных и космических изделий» Московского авиационного института. e-mail: shklyarchuk@list.ru
Application of the one-step Galerkin method for solving a system of ordinary differential equations with initial conditions
© S.V. Russikikh, F.N. Shklyarchuk Moscow Aviation Institute (National Research University), Moscow, 125993, Russia
A nonlinear oscillatory system described by ordinary differential equations with variable coefficients is considered. It is assumed that in the time interval under consideration, the solution of the system is sufficiently smooth - without discontinuities, collisions and bifurcations. From an inhomogeneous system of equations, terms that depend linearly on coordinates, velocities and accelerations and terms that depend non-linearly on these variables are explicitly distinguished. A new approach is proposed for numerical solution by the step method of the initial problem described by such a system of ordinary differential equations of the second order. At the integration step, unknown functions are represented as a sum of functions satisfying the initial conditions: a linear Euler solution and several given correcting functions in the form of polynomials of the second and higher degrees with unknown coefficients. The differential equations at the step are satisfied approximately in the sense of a weak solution by the Galerkin method on a system of corrective functions. Algebraic equations with nonlinear terms are obtained, which are solved by iteration, starting in the first approximation with a linear solution. The resulting solution at the end of this step is used as the initial conditions for the next step. As an example, we consider one homogeneous second-order differential equation without the first derivative with strong cubic nonlinearity in coordinate (at maximum amplitude, the nonlinear force is twice the linear force). This equation has an exact periodic solution in the form of an integral of the energy of a conservative system, which is used to estimate the accuracy of numerical solutions obtained by Galerkin, Runge-Kutta and Adams methods of the second order, as well as by Radau5 and BDF methods at various time intervals (up to 8000 periods of free oscillations of the system) using various constant integration steps (from 0.0025 fractions of a period). At the same time, in the Galerkin method, four identical correction functions were used at each step in the form of polynomials from the second to the fifth degree. It is shown that for large time intervals of calculations, the Galerkin method has a higher accuracy compared to other numerical methods considered. Therefore, it can be used for the numerical solution of nonlinear problems in which it is required to solve them over long time intervals; for example, when calculating steady-state limit cycles of nonlinear oscillations and chaotic nonlinear oscillations with strange attractors.
Keywords: ordinary differential equations, nonlinear systems, initial problem, numerical solutions, one-step Galerkin method
REFERENCES
[1] Demidovich B.P., Maron I.A., Shuvalova E.Z. Chislennye metody analiza [Numerical methods of analysis]. Moscow, Nauka Publ., 1967, 368 p.
[2] Bakhvalov I.V., Zhidkov N.P., Kobelkov G.M. Chislennye metody [Numerical methods]. Moscow, Laboratory of Basic Knowledge Publ., 2000, 630 p.
[3] Collatz L. Chislennye metody resheniya differencial'nyh uravnenij [Numerical methods for solving differential equations]. Moscow, Foreign Literature Publ., 1953, 460 p.
[4] Hairer E., Nersett S., Wanner G. Reshenie obyknovennyh differencial'nyh uravnenij. Nezhestkie zadachi [Solve ordinary differential equations. Nonrigid problems]. Moscow, Mir Publ., 1990, 512 p.
[5] Hairer E., Wanner G. Reshenie obyknovennyh differencial'nyh uravnenij. Zhestkie i differencial'no-algebraicheskie zadachi [Solve ordinary differential equations. Rigid and differential-algebraic problems]. Moscow, Mir Publ., 1999, 685 p.
[6] Skvortsov L.M. Explicit multistep method for the numerical solution of stiff differential equations. Computational Mathematics and Mathematical Physics, 2007, vol. 47, no. 6, pp. 915-923.
[7] Latypov A.F., Nikulichev Yu.V. Numerical methods based on multipoint Hermite interpolating polynomials for solving the cauchy problem for stiff systems of ordinary differential equations. Computational Mathematics and Mathematical Physics, 2007, vol. 47, no. 2, pp. 227-237.
[8] Bulatov M.V., Tygliyan A.V., Filippov S.S. A class of one-step one-stage methods for stiff systems of ordinary differential equations. Computational Mathematics and Mathematical Physics, 2011, vol. 51, no. 7, pp. 1167-1180.
[9] Skvortsov L.M. A fifth order implicit method for the numerical solution of differential-algebraic equations. Computational Mathematics and Mathematical Physics, 2015, vol. 55, no. 6, pp. 962-968.
[10] Belov A.A., Kalitkin N.N. Curvature-based grid step selection for stiff Cauchy problems. Mathematical Models and Computer Simulations, 2017, vol. 9, no. 3, pp. 305-317.
[11] Volkov-Bogorodskii D.B., Danilin A.N., Kuznetsov E.B., Shalashilin V.I. Implicit methods for integration of initial value problems for parameterized systems of second-order ordinary differential equations. Computational Mathematics and Mathematical Physics, 2003, vol. 43, no. 11, pp. 1620-1631.
[12] Bulatov M.V., Gorbunov V.K., Martynenko Yu.V., Kong N.D. Variational approaches to numerical solution of differential algebraic equations. Computational Technologies, 2010, vol. 15, no. 5, pp. 3-13.
[13] Chistyakov V.F., Chistyakova E.V. Application of the least squares method to solving linear differential-algebraic equations. Numerical Analysis and Applications, 2013, vol. 6, no. 1, pp. 77-90.
[14] Zaletkin S.F. Numerical integration of ordinary differential equations using orthogonal expansions. Matematicheskoe modelirovanie [Mathematical modeling], 2010, vol. 22, no. 1, pp. 69-85.
[15] Arushanyan O.B., Zaletkin S.F. On some analytic method for approximate solution of systems of second order ordinary differential equations. Moscow University Mathematics Bulletin, 2019, vol. 74, no. 3, pp. 127-130.
[16] Aulchenko S.M., Latypov A.F., Nikulichev Yu.V. Metod chislennogo integriro-vaniya sistem obyknovennyh differencial'nyh uravnenij s ispol'zovaniem inter-polyacionnyh polinomov Ermita [Method of numerical integration of systems of ordinary differential equations using Hermite interpolation polynomials]. Zhur-nal vychislitel'noj matematiki i matematicheskoj fiziki [Journal of Computational Mathematics and Mathematical Physics], 1998, vol. 38, no.10, pp.1665-1670.
[17] Latypov A.F., Popik O.V. Computational method for solving of the cauchy problem for stiff systems of ordinary differential equations based on multilink interpolated hermite polynomials. Computational Technologies, 2011, vol. 16, no. 2, pp. 78-85.
[18] Ershov N.F., Shakhverdi G.G. Metod konechnyh elementov v zadachah gidrodinamiki i gidrouprugosti [Finite element method in problems of hydrodynamics and hydroelasticity]. Leningrad, Shipbuilding Publ.. 1984, 240 p.
[19] Ioriatti M., Dumbser M. Semi-implicit staggered discontinuous Galerkin schemes for axially symmetric viscous compressible flows in elastic tube. Computers and Fluids, 2018, vol. 167, pp. 166-179.
[20] Kulikov G.Y., Khrustaleva E.Y. Automatic step size and order control in one-step collocation methods with higher derivatives. Computational Mathematics and Mathematical Physics, 2010, vol. 50, no. 6, pp. 1006-1023.
[21] Weiner R., Kulikov G.Y. Efficient error control in numerical integration of ordinary differential equations and optimal interpolating variable-stepsize peer methods. Computational Mathematics and Mathematical Physics, 2014, vol. 54, no 4, pp. 591-607.
[22] Bazhenov V.G., Chekmarev D.T. Reshenie zadach nestacionarnoj dinamiki plastin i obolochek variacionno-raznostnym metodom [Solving problems of unsteady dynamics of plates and shells by the variational-difference method]. Nizhny Novgorod, UNN Publ., 2000, 118 p.
[23] Belytschko T. A survey of numerical methods and computer programs for dynamic structural analysis. Nuclear Engineering and Design, 1976, vol. 37, I ss. 1, pp. 23-34.
[24] Houbolt J.C. A recurrence matrix solution for the dynamic response of elastic aircraft. Journal of Aeronautical Sciences, 1950, vol. 17, pp. 540-550.
[25] Newmark N.M. A method of computation for structural dynamics. ASCE Journal of the Engineering Mechanics Division, 1959, vol. 85, pp. 67-94.
[26] Bathe K.J., Wilson E.L. Numerical methods in finite element analysis. New York, Prentice-Hall, 1976, 544 p.
[27] Krieg R.D., Key S.W. Transient shell response by numerical time integration. International Journal for Numerical Methods in Engineering, 1973, vol. 7, no. 3, pp. 273-286.
Russkikh S.V., Dr. Sci. (Phys. — Math.), Associate Professor, Associate Professor of the Department of Design and Strength of Aeronautical Missile and Space Products, Moscow Aviation Institute (National Research University). e-mail: sergey.russkih@rambler.ru
Shklyarchuk F.N., Dr. Sci. (Eng.), Professor, Professor of the Departmen of Design and Strength of Aeronautical Missile and Space Products, Moscow Aviation Institute (National Research University), Moscow Aviation Institute (National Research University). e-mail: shklyarchuk@list.ru
