Явный метод Эйлера построения последовательности Текст научной статьи по специальности «Математика»

УДК 517.537

ЯВНЫЙ МЕТОД ЭЙЛЕРА ПОСТРОЕНИЯ ПОСЛЕДОВАТЕЛЬНОСТИ

И. Н. Зубов, А. В. Зубов, А. Ф. Зубова,

Н. Н. Учватова

В данной статье изложен метод построения ортогональной матрицы полярного разложения матрицы, использующий аппарат матричных функций Ляпунова. Эта задача построения полярного разложения невырожденной матрицы, т. е. представления ее в виде произведения симметричной положительно определенной матрицы на ортогональную, является весьма важной в приложениях. Такое представление всегда возможно и однозначно.

Рассмотрим линейную систему

X = А ^) X, (1)

где X = {Х1, хп) , А (Ь) — кососимметрическая (п х п) матрица. Условия теоремы существования и единственности мы будем считать выполненными. Решение задачи Коши X = Х0 при Ь = 0 дается формулой

X (Ь) = 3 (Ь) Хо, (2)

где 3 (Ь) — фундаментальная матрица, удовлетворяющая соотношениям

3 = А (Ь)3 3 (°) = Е.

Известно [4], что 3(Ь) является унитарной матрицей, т. е. выполняется соотношение

х (Ь) 3* (Ь) = Е. (3)

Предполагаемая модификация численных методов заключается в таком изменении, что соотношение (3) будет выполнено для всей последовательности приближенных значений 3k фундаментальной матрицы в узлах tk = kh, причем это изменение производится в пределах локальной точности метода, что, очевидно, не ухудшает его [1].

Рассмотрим явный метод Эйлера построения последовательности 3k

3к+1 =3k + ЬА {Ч) 3k. (4)

Поскольку точность этого метода является величиной порядка к2, соотношение

3к+1 = 3к + ЬАк3к + 0 (^2 ), где элементы матрицы О (Ь2) — суть величины порядка Ь2 не «хуже» метода (4). До-

пустим, что соотношение (3) на п-м шаге выполнялось. Тогда имеем

Е„ХП+1 = Е - к2А„2. (5)

Введем этап коррекции в метод (4): вместо матрицы 3п+1 в качестве очередного приближенного значения мы будем брать унитарную матрицу ип+1 полярного разложения матрицы 3п+1 [1; 7].

Пусть хп+1 = ^+1^+1 — полярное разложение матрицы 3п+1, где ^п+1 — симметрическая положительно определенная матрица; ип+1 — унитарная матрица. Из соотношения (5) следует:

^2+1 = Е - к2А2.

Отсюда имеем: Еп+1 = Е - (1/2) Ь2А + О (Ь4 ).

Оценим величину ХП+1 - и„+1. Из последнего соотношения видно, что

3 п+1 — ип+1 = — ^ к Апип+1 — 0 (к )

является величиной порядка к2, т. е. находится в пределах точности метода (4). Обозначим и (3) операцию нахождения унитарной матрицы и полярного разложения матрицы 3, т. е. ип+1 = и (3п+1). Таким образом, мы получили двухэтапный вычислительный алгоритм модифицированного метода Эйлера:

1=1 п+1 = Х п + кАпХп,

Х п+1 = и (1=1 п+1) .

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

© Зубов И. Н., Зубов А. В., Зубова А. Ф., Учватова Н. Н., 2012

Хп+1 Р (Хп-д, Хп-д+1, ■■■, Хп, ■■■, Хп+5)

является методом к-го порядка, т. е. имеем соотношение [5]:

X; = Х (]Щ + О (Ък ).

Тогда из соотношения (3) следует:

X; = Х*+1 = ^ = Е + О (^) ■

Следовательно, Г+1 = £ + О (йк), и операция и ( х у+1) не выводит нас из пределов точности метода.

Итак, вводя этап коррекции

X у+1 = и (х у+1), (6)

мы модифицируем исходный численный метод таким образом, что интеграл (3) сохраняется, т. е. приближенные значения удовлетворяют соотношению (3) [2].

Покажем теперь применение модифицированных методов при интегрировании неоднородных систем. Рассмотрим систему

X = А ^) X + Г X).

По формуле Коши имеем

X (I) = X (I) X-1 (0) Хо + |X (I) X-1 (т) ¥йх. (7) о

Используя формулу прямоугольников, можно получить явный алгоритм

Хк+1 = Xk+iXkXk + Р {гк Хк){Ч+1 - Ч)

и неявный алгоритм

Хк+1 = Xk+1X*kXk + Р (4+1 Хк+1 ){tk+1’tk ), Хк+1 =Хк+1ХкХк + (Р (ЬЬ Хк ) +

+ Г ^+1, Хк+1)) tk^2tk-i.

Можно выписать также предсказываю-ще-исправляющий алгоритм

Хк+1 = ХкХк + Р (^к, Хк )^к+1 - Ь ) ,

Хк+1 =Хк+1ХкХк + Р (¿к+1, Хк+1)( 4+1 - ^к ) ■

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

Рассмотрим матричное дифференциальное уравнение

X = 1 (-X + X *-1). (8)

Приступим теперь к построению численного алгоритма построения матрицы U. Возьмем за основу какой-нибудь одношаговый метод, например метод типа Рунге — Кутты. Погрешность e нахождения ортогональной матрицы U, определяемая выражением e = ||х - U||, где X — приближенное решение задачи, складывается из погрешности 8j — аппроксимации значения U решением дифференциального уравнения (8)

£! = ||х (T, A) - U ,

и S2 — полной погрешности метода численного интегрирования системы (8):

е2 = ||х (T, A) - X.

Величина Si может быть оценена по формуле (9):

V = (V0 - E) exp (-t) + E, (9)

откуда можно получить оценку для требуемой длины интервала T интегрирования. Имеем

ei < 11AA* - E| exp (-T).

Отсюда находим T > ln(||AA*||.

Величина S2 в свою очередь зависит от погрешности r формулы используемого метода и от ошибки округления правой части. Погрешность округления правой части связана с целесообразностью вычисления обратной матрицы итерационными методами.

Теперь получим простые оценки величины шага интегрирования, числа шагов и точности нахождения правой части. Для этого зададимся требуемой точностью нахождения матрицы и будем предполагать, что погрешности Si и S2 вносят одинаковый вклад в погрешность e : Si = S2 = e/2 [6].

Для простоты рассуждений можно потребовать, чтобы влияние каждой локальной погрешности на очередном шаге интегрирования не превосходило величину e 2/N + 1, где N — число шагов. Тогда имеем: e2

r + S < —,

N

где r — погрешность формулы метода; 8 — погрешность округления правой части. При построении одношаговых методов обычно

выполняется требование r = O (hk+1), где

k > 0 (число k называют порядком метода). Отсюда видно, что для согласованности вычислительного процесса величина 8 должна удовлетворять требованию 8 = О (йк+1).

Поскольку величина интервала интегрирования уже нами оценена, для нахождения шага А можно выписать следующее неравенство:

Иначе говоря, шаг А должен быть выбран таким образом, чтобы выполнялось соотношение

= О (ны), 1 ^ 0.

Справедлива теорема.

Теорема 1. Решение системы матричных дифференциальных уравнений

У = 1 (-У + УУ*У) (10)

с начальным условием У§ = А*-1 сходится при Ь ^ да к значению матрицы и полярного разложения (11), причем справедлива оценка (9).

А = FU. (11)

Численный метод интегрирования системы (8) может быть основан на методе типа Рунге — Кутты при условии согласования порядка метода и шага с требуемой точностью вычислений.

Можно предложить и другой способ построения полярного разложения матрицы X. Заметим, что справедливо соотношение

ЕЕ = Р . Следовательно, задача нахождения полярного разложения сводится к нахождению положительно определенного корня из матрицы ЕЕ . Матрица и определится по формуле и = Р X [3].

Обозначим через М множество всех квадратных корней из матрицы XX*, т. е.

М = {ф|ф2 = XX*}.

Теорема 2. Множество М является инвариантным, асимптотически устойчивым в целом для матричной системы дифференциальных уравнений

X = -X + X-1xx*. (12)

Доказательство. Инвариантность М очевидна. Для доказательства асимптотической устойчивости множества в целом построим матричную функцию Ляпунова V = X - ЕЕ . Ясно, что V обращается в нуль-матрицу лишь при X е М и удовлетворяет дифференциальному уравнению V = -2V, откуда видно, что V ^ 0 при Ь ^ да при любом начальном значении V) = Хо - XX*. Следовательно, X2 ^ ЕЕ* независимо от начального значения Xо, а это означает, что

X ^ М.

Х0еМ„

Теорема доказана.

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Зубов А. В. Динамическая безопасность управляемых систем / А. В. Зубов, Н. В. Зубов.

СПб. : Изд-во НИИ химии СПбГУ, 2009. 172 с.

2. Зубов А. В. Математические методы безопасности управляемых систем и методы анализа неста-

ционарных систем управления / А. В. Зубов, Н. В. Зубов, Н. И. Зубов. СПб. : Мобильность плюс, 2010. 319 с.

3. Зубов А. В. Теория устойчивости и теория квазипериодических систем / А. В. Зубов, Н. В. Зубов, С. А. Стрекопытов. СПб. : АООТ Мобильность плюс, 2010. 206 с.

4. Зубов А. В. Анализ систем управления и равновесных движений / А. В. Зубов, И. В. Зубов,

М. В. Стрекопытова. СПб. : Мобильность плюс, 2011. 347 с.

5. Зубов А. В. Математические методы исследования устойчивости и надежности технических систем / А. В. Зубов. СПб. : ВВМ, 2011. 362 с.

6. Зубов Н. В. Безопасность функционирования технических систем / Н. В. Зубов, А. Ф. Зубова. СПб. : ВВМ, 2009. 343 с.

7. Зубов С. В. Анализ равновесных движений и расчетная устойчивость / С. В. Зубов, М. В. Стрекопытова. СПб. : СПбГУ, 2010. 446 с.

Поступила 01.11.2011.

УДК 004.658

АНАЛИЗ УПРАВЛЯЕМОСТИ НЕЛИНЕЙНЫХ АФФИННЫХ СИСТЕМ УПРАВЛЕНИЯ В СИСТЕМЕ МА^АВ

В. В. Афонин

Рассматриваются вопросы управляемости нелинейных по состоянию и линейных по управлению систем, которые принято называть аффинными системами управления. Анализ локальной управляемости определяется на основе вычисления детерминанта или ранга матрицы управляемости заданной нелинейной аффинной системы. Формирование матрицы управляемости осуществляется с привлечением алгебры Ли. Для решения данной задачи используется система математических вычислений МА^АВ, с помощью которой, в частности, выполняется оценка времени заключения об управляемости систем с последовательным соединением нелинейных звеньев.

Х (і) = [х (і)1, х (і)2 , х (і)„]Т —

вектор состояния.

В кратком виде можно записать аффинную систему в следующем виде:

— = А (х) + В (х) и.

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

© Афонин В. В., 2012

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

Класс нелинейных динамических систем линейных по управлению называется аффинными системами [1—2]. Математическая модель аффинной системы со скалярным управлением представляется в виде

^ А (х (і)) + В (х (і)) м (і),

где

А (х )) = [«1 (х ))> а2 (х (Ь))> ■■■> ап (X (Ь))]Г,

В (X (І)) = [Ъ (X (І)),

Ъ2 (X (і)), ..., Ьп (X(і))]т,