Двусторонние методы интегрирования жестких систем обыкновенных дифференциальных уравнений на основе интеграла Дюамеля Текст научной статьи по специальности «Математика»

Сок 10.24411/2409-5419-2018-10076

двусторонние методы интегрирования жестких систем обыкновенных дифференциальных уравнений на основе интеграла дюамеля

арутюнян

Тигран Робертович1 НЕКРАСОВ

Сергей Александрович2

Сведения об авторах:

1магистр Московского технического университета связи и информатики, г. Москва, Россия, rob57@mail.ru

2д.т.н., профессор кафедры прикладной математики Южно-Российского государственного политехнического университета, г. Новочеркасск, Россия, пекгавоА: novoch@mail.ru

АННОТАЦИЯ

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

КЛЮЧЕВЫЕ СЛОВА: двусторонний метод; интервальный метод; погрешность; устойчивость; сходимость.

Для цитирования: Арутюнян Т.Р., Некрасов С.А. Двусторонние методы интегрирования жестких систем обыкновенных дифференциальных уравнений на основе интеграла Дюамеля // Наукоемкие технологии в космических исследованиях Земли. 2018. Т. 10. № 3. С. 64-73. Сои 10.24411/2409-5419-2018-10076

Введение

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

Рассматриваемой проблеме посвящено множество публикаций, из которых можно отметить [1-8].

Суммарная ширина wj = wXj + w2 j может быть вычислена из уравнений:

w^, = w. + hw. + Ih1, j = 0,1,...

j+i j j J

Решение данного рекуррентного уравнения равно:

Wj = 2h[(1 + h)) -1] * 2h[exp(tj)-1], j = 0,1,...

За границу допустимого интервала [-1,1] интервальное решение также выйдет при tmax « -lnh, h ^ 0. Принципиально, что для методов более высокого порядка точности к оценка остается той же:

t ~-lnh'k = -klnh = OHnh"1), h^0.

max 4 "

1. Факторы, обуславливающие погрешность двусторонних и интервальных методов решения обыкновенных дифференциальных уравнений (ОДУ)

Основной недостаток известных интервальных и двусторонних методов для решения задачи Коши [1-2] связан с тем, что их погрешность (ширина полосы решения) растет экспоненциально в зависимости от длины интервала интегрирования (0, /Д Для любой заданной точности £ предельно допустимая длина интервала интегрирования весьма слабо (логарифмически) зависит от значения шага интегрирования к: tFшж(z,к) = О(1пк), к ^ 0. Т.е., чтобы увеличить длину промежутка интегрирования в 2 раза, условно говоря, требуется уменьшить шаг интегрирования в в2 ~ 7,4 раз независимо от порядка точности метода. Рассмотрим следующий характерный пример.

Пример 1. Требуется решить задачу Коши для системы ОДУ

dxl| dt = х2, dx2l dt =-х1, t е(0, tp );

X (0) = 1,х2 (0) = 0.

Ее точное решение х1 =cos(t), х2 =бш(/) принадлежит промежутку [-1, 1].

1. Конечно-разностная схема метода Мура первого порядка имеет вид [1, 2]:

Хи + 1 = X, , + hX2: ; + [-1,1^ / 2, X2, ; + 1 = = X2,] - X+ [-1,1]^/2, ] = 0,1,...

X , Х2.—интервалы, содержащие решение: хк(/.)еХ к= 1,2;_/=1,2,...; к — шаг интегрирования. Ширина интервального решения удовлетворяет уравнениям:

2. Двусторонний метод [2]. Решение отыскивается в виде интервалов

X, = si + [-1,1] s/2) + [-а, а^,(1), I = 1,2;

(1) (2)

где а — некоторая неизвестная константа, si, s> , si ' — эрмитовы сплайны третьей степени. Сплайны , = 1,2, аппроксимируют решение х,, , = 1,2, задачи Коши. Они строятся по данным приближенного решения задачи каким-либо численным методом, например Рунге-Кутта. Сплайны si(1), si<2),I = 1,2 , аппроксимируют решения и,V,i = 1,2, двух вспомогательных задач Коши:

du / dt = Wu + w, t е(0, tF ), и (0) = 0, сЬ / dt = Wv, t е(0, tF), V (0 ) = 0,

где вектор V имеет единичные компоненты, а матрица Ж состоит из элементов

=д% / дх, г = 1,2;Жу = |д£ / дх.|г Ф ], г, ] = 1,2,

где /.,1 = 1,2 — функции правых частей задачи Коши.

Для обеспечения двусторонности оценок Х.,,= 1,2, значение параметра а выбирается по формуле:

а = тах(Ф . (^ / (t); t е (0, tF), г = 1,2).

Функции Ф,(/) и определяются при помощи следующих соотношений:

Ф; (t) = | Ф; (s,t)|-ds^ldt + g^(t) + gias^(t)

4,j+1 "1,j 2,j 1 " ' "2,j+1 = W2,j + hwi j + h2, j = 0,1,...

Ч, (t) = ds?4 dt - (t)-gi,2 s2(1) (t),

#

ф,. (5, t) = f¡ ) - ds¡ / dt, g¡ г =дд1 / дх1,g¡ 7 = = \ddfi / дх-|, г ф у,г,] = 1,2,

С учетом исходных данных задачи / (х)=х2, /2(х) = -х1. Следовательно,

Жп = W22 = 0,Шп = Ж21 = 1;

V,.(?) = 0,s()(t) = 0, Фг(?) = |Фг(*,?)|,

Х¥1 () = ds1l) / dt - s2<1) (t),¥2 (t) = ds2m / dt - (t),

X = s¡ + [-а,а]si(1),Уг е (0, гр),г = 1,2.

Предположим сплайны • , = 1,2, аппроксимируют решение исходной задачи х,, , = 1,2, с порядком точности к. Решение вспомогательной задачи Коши ^^ находится точно: u(t) = в'- 1,, =1,2. Тогда =Т2(') = 1.

Ширина двустороннего решения равна:

wJ = а (в1 -1), 1 = 0,1,...; а = 0(Нк) , h ^ 0.

Следовательно, для рассмотренного двустороннего метода сколь угодно высокого порядка точности к оценка остается такой же неудовлетворительной:

~ 1п к- = -к1п к = 0(1п к'1), к ^ 0.

Применимость метода для случая жестких систем ОДУ

В работе [8] показано, что двусторонний метод [2], в принципе, позволяет решать некоторые жесткие задачи Коши для ОДУ. Однако, соответствующий класс задач относительно узок. Например, в статье [8] приводится пример жесткой системы с матрицей:

А =

-20000 0 0 20000 -1 1

0

1 -1

Однако достаточно изменения значения одного элемента матрицы а23с 1 на -2, как двусторонний метод перестает быть устойчивым, хотя система ОДУ остается асимптотически устойчивой. Отмеченные свойства определяются собственными значениями матриц A и Ж, которые соответственно равны:

\ = -20000, Я2>3 = -1 ± ¡Л, Re\2,з < 0;

Я, =-20000, Я2 =72 -1, Яз =-72 -1, Re Я2 > 0.

3. Традиционные численные методы интегрирования ОДУ

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

Погрешность традиционных численных методов для отмеченных моделей относительно слабо возрастает по мере увеличения интервала интегрирования. Для иллюстрации рассмотрим решение задачи Коши простейшим традиционным методом наименьшего порядка точности— явным методом Эйлера. Его разностная схема имеет вид:

Х1, j+1 _ Х1, j + ^Х2, j, Х2, j +1 _ Х2, j ЬХ1, j ' ] = 04, • • • .

Формула для сеточного решения имеет вид:

= х1,; + 'хи, I2 = -1, 7 = 0,1,...;

= ( -)) 2о, 7 = °4,---;

откуда находится оценка погрешности численного решения:

\в"' - (1 - гй) | < '¥йе'гН <е, у = 0,1,...; и допустимая длина

интервала интегрирования: ^ ю—= О(ПГ1),h ^ 0.

h

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

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

В теории интервальных методов решения ОДУ существует также проблема, носящая название «эффект

раскрутки» («эффект упаковывания», эффект Мура), широко обсуждаемая в литературе [1,2] и имеющая некоторое значение для описываемых ниже методов. Эффект Мура определяется как явление расширения интервального решения динамической системы вне зависимости от применяемого метода (в рамках интервального подхода). Этот эффект обычно иллюстрируется на примере задачи:

dx1 / dt = х2, dx2 / dt = - х1, t > 0,

X (0)е1 + [-е,е],х2 (0)е[-е,е], е> 0.

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

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

2. Интервальный метод для решения

жестких систем ОДУ, устойчивый на всей числовой оси

2.1. Постановка задачи

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

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

Без ограничения общности будем рассматривать задачу Коши для автономной системы ОДУ

dx / dt = f (x),t > 0;x(0) = x0

(1)

где x и f- вектор-функции x = (xp ..., xm)T, f= f ..., fm)T , x0 = (x01, ..., x0m)T, значения начальных условий x0 заданы с погрешностью 5, f (x) = Ax + r(x), A = [aij \ j=l m — матрица с действительными коэффициентами, все собственные значения {Xk }k=1 m, которой имеют отрицательные действительные части, строго меньшие величины — Х0 < 0; r(x) — нелинейная часть вектор-функции fx), обладающая свойством max{|| r(х)||,|| rx(x)|J|x ||} < C01| x ||2, Vx: || x ||<y, где || x || = max{| xl\, xm |}, rx (x) — матрица Якоби [dr / dxj];>j=1, m, || - f. —согласованная с INI матричная норма: | A |1 = max{| ал | + ...+ | ain | ;i = 1, ...,m}, C0, у — некоторые положительные константы.

Будем предполагать, что решение задачи (1) существует, единственно и принадлежит классу C1 (0, ° ), причем x(t) ^ 0 при t ^ да. Следовательно, существует конечный m — мерный интервал А, такой что x(t) е А Vt > 0.

Введем в рассмотрение разбиение полуоси (0, ° ) с шагом h: 0 = t0 < tx < ...< tk < ....

При помощи интеграла Дюамеля задача Коши (1) редуцируется к эквивалентному интегральному уравнению Вольтерра:

x(t) = exp [ A(t-tj_x) ] x(tj_x) + t

+ J exp [A(t- t)](x(T))dt ,

(2)

j-1

I > ] = 1,2,... .

Использование этой формулы может приводить к вычислительным методам двух типов, рассматриваемых в последующих частях статьи.

2.2. Многошаговый интервальный метод на основе решения интегрального уравнения Вольтерра

Это уравнение получается из (2) при ] = 1:

t

х() = ехр(At)х0 +1 ехр - т)] г(х(т))^т , (3) 0

Интервальные оценки решения X. в узлах интегрирования t = tj находятся в соответствии с правилами интегрирования интервальнозначных функций [1] по формулам:

Xj = exp(Atj)X0 + £exp[A(L-Tk)]R(X* )hk

(4)

k=i

где . = 1,2, ...; Х0 — интервальная оценка начального условия х0, Хк" - интервальная оценка решения х(') на отрезке Тк = [к-1, tk ], \ = tк - tk-1, Я(X) - интервальное расширение функции г(х), для которого будем предполагать выполненным условие: Х))|| < С||Х||-||и>(Х)||, где

С — некоторое положительное число. Это предположение согласуется со свойствами функции г(х).

Если известна интервальная оценка X. 1 решения задачи х(') при ' = ' то, согласно (2), интервальная оценка х(') на отрезке Т. (т. е. значение X*) при достаточно малом шаге интегрирования может быть вычислена посредством решения интервального уравнения:

X/ = [елн< X- + И^' R(X;)] п А,И = [0,^],

1 ^ 1 -1 1 1 = 1,2,...

■¡ь (5)

Утверждение 1. Оценки, получаемые при помощи метода (4), (5) обладают свойством двусторонности.

Доказательство относительно элементарно и осуществляется на основе стандартных технических приемов [1,2].

Утверждение 2. Интервальный метод (4), (5) при И ^ 0 сходится и асимптотически устойчив.

ДОКАЗАТЕЛЬСТВО. Оценим сверху значения ширины интервального решения м>(Ху), м>(XУ), у = 1,2,.... Из соотношений (4), (5) с учетом формул для ширины произведения интервалов [1]: w( АВ) < \А\ w(B) + |В| w( А), АН^) < 21A\hj,находим

ЦX,) < \гм> ЩX0)|| +

+± [| „ (^)х* )|| + || X* ))

(6)

Х/)|| < \е >\ -||м>(Х]_1)| +

,, ........................(7)

+ иеШ<) Лх.Л + 2к]\еШЛ\ -ад*),

] = 1,2,...

Так как действительные части всех собственных значений матрицы A отрицательны и решение задачи Коши (1) ограниченно и стремится к нулю при t ^ да, то обоснованы оценки:

\\еА'\\ < С1е~и, || х(1 )||< СеХ1,

(8)

где X, С — некоторые положительные константы. Используемые ниже коэффициенты С2, С3, ... также обозначают некоторые положительные числа.

Идея доказательства основана на стандартном техническом приеме: вначале предположим, что для интервального решения (4), (5) выполнены оценки, аналогичные (8):

11X11 ^ Св'

(9)

далее на основании соотношений (6)-(9) получим оценки для ширины интервального решения и покажем, что они согласуются со сделанными предположениями.

Так как и>(Х0) <8,

R

(X**))

< с X*

к)

< С, X*

\\еАНк\\ < С

С учетом (5):

\\х*\\ <|\еАНк\\(\\Хк-Л + hkR(Хк*))< <|\еАНк\ |(|| Хк + hkС\\Хк | )<

< С2 (||хкЛ + иДЦхкII2) < С21|Хк Л + С3\ ||Х,

2

к II '

следовательно, при достаточно малых кк:

||х;|| < С4||хк-Л < С5^<к, к)

< СК,

С учетом (7) ||Ч X, *)|| < С7 ((X, + hk),

К X. )|| <С/Л' 5 + С8 £

¡, 2е'2Х11 + ¡,

< €,е"> 5 + С, X[¡,,V' + Ь,\„(X,-! + К V' ],

к=1

После подстановки (8)-(9) в (6) и некоторых преобразований находим:

X )| < С10 ехр(-Ц) 1 = 1,2,...

1 -1

8 + h + Хк I, (10)

где к — наибольший шаг интегрирования.

Из соотношений (10) после ряда технических преобразований может быть получена оценка, которая согласуется со сделанными предположениями (9) при 5, к ^ 0:

)|| < С4ехр Щ) (8 + h), ] = 1,2,... (11)

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

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

к=1

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

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

2.3. Одношаговые интервальные методы

Его простейший вариант имеет вид:

гДе aiMi е {0,1},i = 1,...,m -1,

X. = e >X._l + h.eAH'R(X '), j = 1,2,

(12)

Интервальные оценки X * решения на отрезках [^_1, tj ], у = 1,2,... , в данном методе вычисляются также при помощи соотношений (5). Из (12) находится следующая оценка ширины интервального решения:

\\w(X.)|| < I\eM' II • ||w(X-J| + к-\\w(eAHR(X'))\\ | j = 1,2,

Поскольку ||е 71| может быть больше 1, то данный вариант метода, вообще говоря, не является асимптотически устойчивым, что подтвердили вычислительные эксперименты для традиционных модельных задач.

Отмеченный недостаток метода удается эффективно устранить при помощи специальной замены переменных.

В дальнейшем, без ограничения общности шаг интегрирования считается постоянным.

Пусть D = D(h) — преобразующая матрица, соответствующая жордановой форме Л матрицы B = eAh: Л=D"1BD. Согласно известной теории матриц существует преобразующая матрица D такая, что жорданова форма Л имеет вид:

Л =

exp (\h) 0

0 exp (X 2 h) a23^

n = q -max{ehR*k;k = 1,...,m} =

= e~x°h - ehRXi > 0 ,

|< q = e-x'h < 1.

Матрицу Л можно и, как правило, предпочтительней выбирать диагональной, полагая ац+1 = 0,1 = 1, —, т -1, при этом матрица D состоит из собственных векторов матрицы A (каждому /-тому столбцу матрицы D соответствует собственный вектор, отвечающий .-тому собственному значению).

Осуществим замену зависимых переменных по пра-вилу:}() = D'lx(í), тогда из (2) следует, что у(/=+ , / = 1,2, ...,

']

8] = 1 ехР[А('] -т)]г(х(т))т.

']-1

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

Л () = У0к, 1 + Рк ()ехри<?к ()

к = 1, ...,т; 1 = 1, 2 ...;i2 =-1,

о

где значения ук ] находятся при помощи соотношений:

уЬ=X Л »У—+к=1,..., т; з=1,2,.•..

1=1

где Б- = №{г1)к1, ]к1=Кт1 ,шld(Glj)—середина интервала О.

^ = Ь^'И ) ^ = (,, .

Тогда

Р* ) ехР0фк (tj)] = X Л *; Р; (^-1 )ехР[г'Ф; (tj-l)] +

;=1

(g;,, - mid(Gi/)),к = 1,...,т;] = 1,2,....

Учитывая, что ^ j - mid(Gг ;)| < 1/2 ;), где м>(О1 ) — ширина интервала О, , находим рекуррентные соотношения для вычисления оценок сверху к значениям

РД):

РК! = Х|Лк I |РО-1 +1/21^ ),

I=1

к = 1,..., т;! = 1,2,....

(14)

Интервальные оценки X к решению задачи Коши х(') могут быть вычислены по формуле:

т

Хк1 =£ Re ((А;) + [-1,1]| ри,

(15)

к = 1,..., т;; = 1,2,..

Начальные значения оценок вычисляются при помощи соотношений:

У0 = П-^(Х0), ры = 1/2ЦX,,0), (16) к = 1,..., т;х0 е Х0.

Утверждение 3. При сделанных в п. 2.1 предположениях относительно свойств задачи Коши (1) интервальное решение (14)-(16), также как точное, стремится к нулю при да, а его погрешность удовлетворяет оценке Цр,. || = 0(8 + И) при 5, к ^ 0, где 5 — погрешность задания начального условия задачи Коши.

ДОКАЗАТЕЛЬСТВО. Оценим по норме ширину интервалов 01 .:

)|| < Сп (21X)|| + h||w(R(X]))||) < < С12 (Х]\|2 + ^|х]||||цX])||)< С,3( + X4),

Так как w(X.) = р., то из (14) следует:

|Р || < д ||р,+1/21|\wiGj )|| < д ||р+

+ с14( ( + р;) = (д + С14 к )) + С14 к2, } =

где при достаточно малом к: д1 = д + С14к< 1, следовательно,

||р,|| ^ С15^ + Cl6k, У =12 • • •, что и доказывает утверждение.

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

Повышение порядка точности данных методов возможно при помощи разложения в (2) вектор-функции г(х(т)) в ряд Тейлора по степеням (т - ') с интервальной оценкой остаточного члена.

Пример 2. Рассмотрим систему ОДУ, описывающую переходной процесс в нелинейной электрической цепи:

Ldi / dt + ге11 + ис = 0, 1 = dq / dt, ис = q / С^),^ > 0;

(17)

I (0)= 0, ыс (0)= ыс0, С()= С0/ (1 + kq2 ). (18)

В результате преобразований система уравнений (17)-(18) приводится к задаче Коши для дифференциального уравнения Дюффинга, описывающего нелинейные колебания [10]:

d2г / dт2 + а dz / dт + г + г3 = 0, т> 0; г(0)= q0к1'2, dz / dт(0)= 0,

д () = г(г (Щ )"1/2^-1'2,

а = Ге1 (С0/1 )1/2, ис0 = д»/ с(д»).

В нормальной форме постановка задачи Коши имеет вид:

dxl/dх = х2, dx2 /dт = -ах2 - х1 - х13, т > 0; X; (0)= д0к1/2, х2 (0) = 0.

где г = х1, dzldт = х2. Задача из примера 1 является частным случаем данной задачи Коши, так как описывает идеальные периодические колебания. Таким образом,

А =

0 1

-1 -а

г (х) =

0

R( х) =

0

- X3

\Х.

Собственные значения матрицы А равны:

Х1 = [-а - (а2 - 4)1/2]/ 2, = [-а + (а2 - 4)1/2]/ 2.

Матрицу ех^АО находим по формуле Сильвестра [10]:

ехр (At) = [ ехр(Х^)(А -X 2 Е) -- ехр(Х21)(А -Х1 Е)]/(Х1 -Х2).

Следовательно,

ехр(-ак / 2) В = ехр( ЛИ) = ——-1 х

Р

а$Ыкр / 2) + рсЬ(кр / 2) 2вЬ(кр / 2)

-2вЬ(кр / 2) -а$Ъ(кр /2) + рсЬ(кр / 2)

Р

= (а2 - 4)

Преобразующая матрица, соответствующая матрице В = ехр(А^), может быть выбрана в виде:

D =

1 1

[exp(A^h) -bn]/bi2 [exp(^2h) -bn]/bn

Рассмотрим случай а = 531/2, соответствующий апериодическому разряду конденсатора. Собственные значения матрицы А равны: А,= -0,140054945, = -7,140054945. Число жесткости равно А,2 / А,1 ~ 51.

Выбор начального приближения для интервальных оценок решения может быть осуществлен при помощи первого интеграла системы ОДУ, из которого следует, что: ^ (т)| < ^ (0 ), < / ат\< ^ (0 ) [1 +1/2 z2 (0 )]1/2.

Графики интервального решения при различных значениях шага интегрирования и виде начальных условий представлены на рис. 1-3 и качественно иллюстрируют

X(t)>

0.25"

0.20 "

0.15" 0.10" - \\xi(t)

0.05"

0 6 30 t

-0.01" -0.02" / /Ш)

-0.03"

-0.04"

-0.05"

Рис. 1. Интервальное решение при гр = 30; к = 0,3

степень влияния величины шага на точность аппроксима- x(t)'

ции и устойчивость численного решения. 0.25"

При г = р к = 0,3: 0.20" 0.15

Х1 = [0.00008; 0.00410], Ц X) = 0.00418; 0.10" - \ X1(t)

Х2 = [-0.00058; -0.00001], Ц Х2) = 0.00059. 0.05

0 6 30 t

При г = гр; к = 0,15: -0.01"

-0.02"

Х1 = [0.00286; 0.00373], ЦХ) = 0.00088; -0.03" -0.04"

Х2 = [-0.00052; - 0.00040], и<Х2) = 0.00012 -0.05"

Как следует из результатов теоретического анализа и вычислений решение рассматриваемой задачи Коши содержит быстро убывающую компоненту, соответствующую собственному значению А,2. Начиная с т ~ 0,5 решение системы почти полностью определяется медленно убывающей составляющей.

Ограничение на величину шага интегрирования по критерию устойчивости при использовании методов Рун-ге-Кутта определяется неравенством [10] к < 2/ \Х21« 0.28 . Для интервальных методов [1], в силу их известных свойств, имеют место аналогичные или более жесткие ограничения на величину шага.

Как показывают, представленные на рис. 1 результаты, устойчивость численного решения сохраняется и при Н > 2/2|.

По характеристикам вычислительной эффективности для случая тестовой задачи (17) предложенный метод существенно превосходит интервальные и двусторонние методы [1-2]. Например, при помощи интервального метода [1] того же порядка точности решение на рассматриваемом

Рис. 2. Интервальное решение при гр = 30; к = 0,15

x(t)' 0.30 " \

0.24 " к \ \\

0.18" 0.12" v X1(t)

0.06"

0 -0.01" -0.02" 6 / / ^X2(t) i--

-0.03" / /

-0.04"

-0.05"

Рис. 3. Двустороннее решение при гр = 30; к = 0,15 и интервальном начальном условии х1(0)е [0,25; 0,3], х2(0)=0

71

промежутке не удается получить даже при значениях шага интегрирования много меньших 10-6.

Литература

1. Калмыков С. А., Шокин Ю. И., Юлдашев З. Х. Методы интервального анализа. Новосибирск: Наука, 1986. 224 с.

2. Добронец Б. С., Шайдуров В. В. Двусторонние методы. Новосибирск: Наука, 1990. 208 с.

3. Рогалев А. Н. Границы множеств решений систем обыкновенных дифференциальных уравнений с интервальными начальными данными // Вычислительные технологии. 2004. Т. 9. № 1. С. 86-94.

4. Рогалев А. Н. Вычисление гарантированных границ множеств достижимости управляемых систем // Автометрия. 2011. Т. 47. № 3. С. 100-112.

5. Рогалев А. Н. Вопросы реализации гарантированных методов включения выживающих траекторий управляемых систем // Вестник Сибирского государственного аэрокосмического университета им. академика М. Ф. Ре-шетнева. 2011. № 2 (35). С. 54-58.

6. Moore R. E, Kearfott R. B., Cloud M. J. Introduction to interval analysis. Philadelphia: SIAM, 2009. 223 p.

7. Добронец Б. С. Интервальная математика. Красноярск: КГУ, 2004. 216 с.

8. Dobronec B. S. Two-Sided Solution of ODES via a posteriori error estimates // Journal of Computational and Applied Mathematics. 1988. Vol. 23. No. 1. Pp. 53-61.

9. ВолковЕ. А. Численные методы. М.: Наука, 1987.

248 с.

10. Корн Г., Корн Т. Справочник по математике. М.: Наука, 1978. 832 с.

11. Nekrasov S. A. A bilateral method of solving Cauchy problems // USSR Computational Mathematics and Mathematical Physics. 1986. Vol. 26. No. 3. Pp. 83-86.

12. Nekrasov S. A. The construction of two-sided approximations to the solution of a cauchy problem // USSR Computational Mathematics and Mathematical Physics. 1988. Vol. 28. No. 3. Pp. 23-29.

13. Nekrasov S. A. Bilateral methods for the numerical integration of initial- and boundary-value problems // Computational Mathematics and Mathematical Physics. 1995. Vol. 35. No. 10. Pp. 1189-1202.

14. Nekrasov S.A. Efficient two-sided methods for the Cauchy problem in the case of large integration intervals // Differential Equations. 2003. Vol. 39. No. 7. Pp. 1023-1027.

15. Nekrasov S. A. A numerical method for solving dynamical systems with lumped parameters which accounts for an input data error // Journal of Applied and Industrial Mathematics. 2016. Vol. 10. Issue 4. Pp. 528-537.

BILATERAL METHODS OF INTEGRATING STIFF SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS ON THE BASIS OF DUHAMEL INTEGRAL

TIGRAN R. ARUTYUNYAN,

Moscow, Russia, rob57@mail.ru

KEYwORDS: bilateral method; interval method; error; stability; convergence.

SERGEY A. NEKRASOV,

Novocherkassk, Russia, nekrasoff_novoch@mail.ru

ABSTRACT

In work the factors causing an error of bilateral and interval methods of the solution of ordinary differential equations are considered. Examples are given to illustrate the exponential nature of the error growth of the considered methods. Comparison with traditional real methods is carried out. An interval method for solving rigid ordinary differential equations systems that is stable on the whole numerical axis. Bilateral methods with a posteriori error estimation with significantly improved computational qualities for the case of large integration intervals are proposed. Describes an interval method for

numerical integration of initial value problems, which is effective for a practically important class of rigid ordinary differential equations systems and has convergence and stability on an infinite interval. The Cauchy problem for an Autonomous nonlinear ordinary differential equations system is considered without limitation of generality. Values of initial conditions are given with an error. It is assumed that the solution of the initial problem exists only and belongs to class C1. Duhamel integral is used to construct the computational method.The use of this formula leads to two types of computational

methods: a multistage interval method based on the solution of the Volterra integral equation.

Convergence and asymptotic stability of the numerical method with respect to perturbations in initial conditions are proved. Another interval method, derived from the Duhamel equation, is a one-step finite difference method in which the interval estimate of the solution at each step is calculated recursively by the values of the corresponding estimate at the previous step. With equal accuracy, the methods have a number of advantages and disadvantages. The first method is algo-rithmically simple, but the time spent using it increases in proportion to the square of the number of integration steps. The second method is characterized by an order of magnitude smaller total computational cost, linearly dependent on the number of steps of integration as in conventional (i.e. real) methods of integration of odes. However, in General, to ensure the convergence and stability of the second method on the infinite integration interval, a special replacement of variables is required, which complicates the algorithm of calculations. The second method is an interval generalization and development of the so-called semi-analytical method used for the calculation of electrical circuits in the case where the time constants of different parts of the circuit are significantly different. Despite the additional difficulties associated with the calculation of the exponential matrix, this method provides a multiple time gain in practical calculations.

REFERENCES

1. CalmykovS. A., Shokin Y. I., Yuldashev Z. H. Metody interval'nogo analiza [Methods of interval analysis]. Novosibirsk: Science, 1986. 224 p. (In Russian)

2. Dobronets B. S., Shaidurov V. V. Dvustoronnie metody [Bilateral methods]. Novosibirsk: Science, 1990. 208 p. (In Russian)

3. Rogalev.A.N. Bounds of solution sets of differential equations with interval initial data. Computational Technologies. 2004. Vol. 9. No. 1. Pp. 86-94. (In Russian)

4. Rogalev.A.N. Calculation of guaranteed boundaries of reachable sets of controlled systems. Optoelectronics, Instrumentation and Data Processing. 2011. Vol. 47. No. 3. Pp. 287-296.

5. Rogalev.A.N. Implementation of contributed methods of contain-

ment of viable trajectories of controlled systems. Vestnik Sibirskogo gosudarstvennogo ajerokosmicheskogo universiteta im. akademika M.F. Reshetneva [Bulletin of Siberian state aerospace University]. Academician M. F. Reshetnev. 2011. No. 2 (35). Pp. 54-58. (In Russian)

6. Moore R. E., Kearfott R. B., Cloud M. J. Introduction to interval analysis. Philadelphia: SIAM, 2009. 223 p.

7. Dobronets B. S. Interval'naya matematika [Interval mathematics]. Krasnoyarsk: KGU. 2004. 216 p. (In Russian)

8. Dobronec B. S. Two-Sided Solution of ODES via a posteriori error estimates. Journal of Computational and Applied Mathematics. 1988. Vol. 23. No. 1. Pp. 53-61.

9. Volkov E.A. Chislennye metody [Numeric's methods]. Moscow: Nauka, 1987. 248 p. (In Russian)

10. Korn. G., Korn T. Spravochnikpo matematike [Handbook of mathematics]. Moscow: Science, 1978. 832 p. (In Russian)

11. Nekrasov S. A. Two-way method of solving Cauchy. SSSR problems computational Mathematics and mathematical Physics. 1986. Vol. 26. No. 3. Pp. 83-86.

12. Nekrasov S. A. Construction of bilateral approximations to the solution of the Cauchy problem. USSR computational Mathematics and mathematical Physics. 1988. Vol. 28. No. 3. Pp. 23-29.

13. Nekrasov S. A. Bilateral methods of numerical integration of initial-boundary value problems. Computational Mathematics and mathematical Physics. 1995. Vol. 35. No. 10. Pp. 1189-1202.

14. Nekrasov S. A. Effective two-way methods for the Cauchy problem in the case of large integration intervals. Differential equations. 2003. Vol. 39. No.7. Pp. 1023-1027.

15. Nekrasov S. A. A numerical method for solving dynamic systems with lumped parameters, taking into account the error in the input data. Journal of applied and industrial mathematics. 2016. Vol. 10. Issue 4. Pp. 528-537.

INFORMATION ABOUT AUTHORS:

Arutyunyan T. R., Master, Moscow technical University of communications and Informatics;

Nekrasov S. A., PhD, Docent, Professor of applied mathematics, South-Russian state Polytechnic University.

For citation: Harutyunyan T. R., Nekrasov S. A. Bilateral methods of integrating stiff systems of ordinary differential equations on the basis of Duhamel integral. H&ES Research. 2018. Vol. 10. No. 3. Pp. 64-73. doi: 10.24411/2409-5419-2018-10076 (In Russian)