Численное решение систем ОДУ с постоянными коэффициентами

Тип работы:
Реферат
Предмет:
Физико-математические науки


Узнать стоимость

Детальная информация о работе

Выдержка из работы

УДК 519. 622. 2
Численное решение систем ОДУ с постоянными коэффициентами
Д. А. Ахременков, К. П. Ловецкий
Кафедра систем телекоммуникаций Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117 198, Россия
В работе сравниваются несколько подходов к решению систем дифференциальных уравнений с постоянными коэффициентами над полем комплексных чисел, приводятся сравнительные результаты численного и аналитического решения для одной из таких систем.
Ключевые слова: системы дифференциальных уравнений с постоянными коэффициентами, численные методы решения.
1. Методы решения систем дифференциальных уравнений
Рассмотрим систему линейных дифференциальных уравнений с постоянными коэффициентами
? = Ау, у (0) = уо, (1)
где у (х) = (ух (х), у2(х),… уп (х))Т — вектор неизвестных функций, х — независимая переменная- А — (п х п)-матрица постоянных коэффициентов, уо — вектор начальных условий.
В зависимости от характеристик матрицы, А существует несколько подходов к решению системы.
1.1. Все собственные значения различны, собственные вектора
независимы
В случае, если для рассматриваемой системы известны все собственные значения (Ах, …, Ап) матрицы А, они различны и для этих п значений можно выбрать п линейно независимых векторов (ух, Ь2,…, ьп), то, принимая во внимание, что у (х) = V ¦ еХх, получим Ау = XV, откуда следует
/ч т
А = (г& gt-l, г>-2,…, г>-n)diag (АьА2,…, Ап),
что можно представить в виде
Т-1АТ = ^(А1,А2,…, Ап), (2)
где Т — матрица, столбцы которой — собственные вектора матрицы А.
Если применить преобразование у (х) = Тг (ж), = Тг (ж), то первоначальная система дифференциальных уравнений (1) примет вид
= ^ ((Ах), (А2),…, (Ап)) я- х (0) = ?0 = Т-1уо, которую можно легко проинтегрировать:
х (ж) = diag (ехр (Ахж), ехр (2х),…, ехр (Ап ж)) х0.
Авторы выражают благодарность профессору Севастьянову Л. А. за полезные замечания и плодотворное обсуждение статьи.
Окончательно получаем:
у (ж) = Т diag (ехр (Ахж), ехр (А2Ж),…, ехр (А"х)) Т-1уо. (3)
Однако вышеописанная теория хороша, но в ней есть несколько подводных камней:
а) не у всякой матрицы, А существует п линейно независимых собственных векторов, поэтому вышеописанная процедура не применима для таких матриц-
б) даже если такое разложение существует, матрица Т может быть плохо обусловлена, поэтому численная реализация процедуры будет неустойчивой.
1.2. Приведение к диагональному виду
В соответствии с теоремой об аппроксимации [1] для любой квадратной матрицы, А существует матрица Т, такая что:
Т-1АТ =
А1 & amp-12
А2
Ь
1 п
А
0
п
причём ^^ |bi, j| ^ е, где е — произвольное малое положительное число. г, 3
Если? мало по сравнению с машинной точностью, то элементами выше диагонали при численном решении можно пренебречь и решение исходной системы находится так же, как в случае, рассмотренном в пункте 1.1. Однако результат этой теоремы противоречит на первый взгляд тому факту, что произвольная матрица, имеющая кратные собственные значения, не может быть приведена к диагональной форме. Разъяснение сводится к тому, что Т зависит от е. Попытка устремить? к нулю приведёт либо к тому, что матрица Т станет вырожденной, либо к тому, что соответствующий предел перестанет существовать.
1.3. Разложение Шура
Известно [2] (теорема Шура), что для любой комплексной матрицы, А существует такая унитарная матрица Я, что:
(Л1
Я * АЯ =
X X ¦ ¦ ¦ X
А2 X ¦ ¦ ¦ X
СО X
V
/
После приведения матрицы, А к форме Шура и применения преобразования у (ж) = Яу (х), (ж) = Я^'- (ж), получим:
А
П
/ ?1 /Ах Ъ2
(1х
1
Ьх. п-1 Ьхп


'-п- 1, п
1
/ V
(4)
го = г (0) = Я * уо.
(1гп
Последнее уравнение = Хпхп системы (4), можно проинтегрировать независимо от других. В результате получим хп = ехр (Хпх) гп0. Тогда для гп-1 имеем:


Ап-1%П-1 + 1уП^П!
(5)
где хп — известно. Это линейное уравнение (неоднородное, если Ьп-1,п = 0) можно решить методом Эйлера. Возможны два случая:
1. Если Хп-1 = Хп, положим хп-1 = Еехр (Хп-1х) + Еехр (А"ж), подставим это выражение в (5) и сравним коэффициенты. В результате получим Е = Ьп-1,п%по/ (Хп — Хп-1) и Е = гп-1,о — Е-
2. Если Хп-1 = Хп положим гп-1 = (Е + Ех) ехр (Хпж), получим Е = Ьп-1,пгп0 и Е = гп-1,о.
Далее по индукции находим хп-2 и для остальных переменных проводим аналогичные рассуждения.
Для случая XI = Х^ можно получить простую рекурсивную формулу в явном виде:
(х) = X] ехр (^зх)
3 = 1
(6)
Для этого подставим (6) в (4) и, после приравнивания коэффициентов, получим для г = п, п — 1, п — 2, и т. д. [2]
1
Е*к = ^ - Л. Ез^, к = ^ + М + 2,
Ец = гм — Е^ 3=г+1
В случае же равных собственных значений Xi = Х^ простых рекурсивных формул получить не удаётся и в этом случае следует каждый раз повторять процедуру пункта 2).
Окончательно решение исходной системы (1) получается следующим образом
у (х) = - Я г (х).
а
X
г
г
п
п
п
2. Численные эксперименты
На основе разложения Шура разработана программа для численного решения системы (1). Чтобы проверить точность численного решения по этому алгоритму,
расчёты проводились для матрицы А, допускающей аналитическое решение:
/ 0 1 0 0
2 1 -1 -1
0 0 0 1
V -1 -1 2 1)
А =
Эта система имеет следующие собственные векторы и соответствующие им собственные значения:
— 3 — -0. 22 -0. 66 0. 22 0. 66
1 -0.5 -0.5 -0.5 -0. 5
-1 -0. 69 0. 69 0. 14 -0. 14
-1 -0. 39 0. 39 -0. 59 0. 59
Поскольку собственные векторы линейное независимы, общее решение системы будет иметь вид:
(7)
'- Уз = C1v11eXlX + С2^3еА2Х + С3^3еХзХ + C4^3eX4X
У2 = Civ2eXlX + C2veX2X + C3vleX3X + C^veX4X
Уз = Cv3eXlX + C2veX2X + C3 v33eX3X + C4^3eX4X
, У4 = CiveXlX + C2veX2X + C3vieX3X + C4veX4X
Коэффициенты С находятся из начального условия: Yq Подставляя х = 0 в (7) получим:
у0 = Cvl + С2 vi + C3V3 + Ol у0 = Civl + С2 v2 + C3v23 + C4vj у33 = ClV3 + С2 vi + C33VI + CAv3 У0 = ClVf + С2 v4 + + C4vj
{1, 0. 5, -0. 5, -1}.
или в матричной форме: уо = Су, откуда находим вектор коэффициентов С: С = уоУ-1
-1. 705 0
-0. 904 -0. 003
С =
Подставив эти коэффициенты в (7), получим решение системы для заданных начальных данных.
Для той же самой задачи с указанной выше матрицей, А было найдено численное решение для заданных начальных значений. График численного решения приведён на рис. 1.
Решения сравнивались на интервале от 0 до 1, ошибка вычислялась по формуле:
5 = тах щ (ж) — уг (ж), ?=1…4 ®?[0−1] У У Л
где уг (ж) — аналитическое решение в точке ж, а уг (ж) — численное решение в этой же точке. Ошибка составила 8,978Е-13, что является хорошим результатом, учитывая, что аналитическое решение представляет сумму экспоненциальных функций. График ошибки приведён на рис. 2.
Рис. 1. Численное решение системы
1. 20Е-13 1,00 Е-13 3,00 Е-14 6. 00Е-14




д

о о о о о о о о о о о о о о о о о о о о о о о о о
Рис. 2. Отклонение численного от аналитического решения
3. Заключение
Описанный выше метод позволяет решать системы с кратными собственными значениями в комплексной плоскости, однако, хотя алгоритм нахождения такого решения известен, программирование такого решения, само по себе, является нетривиальной задачей.
Литература
1. Беллман Р. Ведение в теорию матриц. — М.: Наука, 1969. — 235 с. [Bellman R. Vedenie v teoriyu matric. — M.: Nauka, 1969. — 235 s. ]
2. Hairer E, Wanner G. Ordinary Differential Equations. — 1993. — 70 p.
UDC 519. 622. 2
Solving Systems of Linear Differential Equations with Constant Coefficients D. A. Ahremenkov, K. P. Lovetskiy
Telecommunication Systems Department Peoples Friendship University of Russia Miklukho-Maklaya str., 6, Moscow, 117 198, Russia
The paper demonstrates the comparison of some different algorithms for solving systems of ordinary differential equations with complex constant coefficients describing the evolution of monochromatic linearly polarized plane electromagnetic waves in a stratified medium.
Key words and phrases: systems of linear differential equations with constant coefficients, numerical methods.

ПоказатьСвернуть
Заполнить форму текущей работой