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



Математические методы. Моделирование. Экспериментальные исследования -►

УДК 517.929

А.В. Зубов, Л.Г. Каляда, А.И. Нечаев, Н.С. Ужегов

МОДИФИКАЦИЯ ЧИСЛЕННЫХ МЕТОДОВ ИНТЕГРИРОВАНИЯ

A.V. Zubov, L.G. Kalyada, A.I. Nechaev, N.S. Ugegov

THE MODYFICATION OF CALCULATING METHODS INTEGRATION

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

система. матрица. соотношение. порядок. точность. значение. множество.

In giving article is supposes the new modification method calculating integration of systems ordinary differential equations is obtains cross-symmetrical matrix coefficients by linear numbers. system. matrix. correlation. order. precision. significance. multitude.

Предложим модификацию численных методов интегрирования систем обыкновенных дифференциальных уравнений, учитывающую физическую природу этих уравнений, т. е. наличие интегралов. Мы будем модифицировать численные методы интегрирования систем обыкновенных дифференциальных уравнений, обладающих кососимметрической матрицей коэффициентов при линейных членах. Модифицированные методы более точны, чем исходные методы Рунге — Кутты, Адамса — Штер-мера и др., и более предпочтительны в задачах интегрирования на больших промежутках. Точность модифицированных методов обусловлена тем, что они сохраняют интеграл однородной системы — свойство ортогональности фундаментальной матрицы.

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

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

X = А(()Х , (1)

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

X () = W(t )Хо, (2)

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

= А(/W(0) = Е.

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

W(t Щ*^) = Е. (3)

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

Wk+1 = Wk + кА^ . (4)

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

Wk+1 = Wk + кА к Wk + 0(к2),

где элементы матрицы 0(к ) суть величины порядка к2, не «хуже» метода (4).

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

^п+Х+1 = Е - к2ЛI.

(5)

Введем этап коррекции в метод (4): вместо матрицы в качестве очередного прибли-

женного значения мы будем брать унитарную матрицу ип+1 полярного разложения матрицы ^,+1. Пусть Wn+1 = Рп+1ип+1 — полярное разложение матрицы ^,+1, где Бп+1 — симметрическая положительно определенная матрица, ип+1 — унитарная матрица. Из соотношения (5) следует

Р„2+1 = Е - к2Л

Отсюда имеем Бп+1 = Е -I 1

п+1

Оценим величину Wn+l - ип+1. соотношения видно, что

к 2Л , + 0(к 2).

Из последнего

Wn+l - и,+1 = - 2 к2Л пип+1 - 0(к4)

является величиной порядка к2 , т. е. находится в пределах точности метода (4).

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

\¥п+1 = Wn + кЛ п Wn

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

W 1 = W 1 W W )

п+1 V п-д' п-д+1' ' "п' ' п+s'

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

WJ = W( ]к) + 0(кк).

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

WJ = WJ = FJ!+l = Е + 0(кк).

Значит, FJ■+1 = Е + 0(кк), и операция U(WJ■+1) не выводит нас из пределов точности метода. Итак, вводя этап коррекции

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

Рассмотрим систему

X = Л(/ )Х + , X).

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

(

X(/) = W(t)W-1(0)Х0 +1 W(t)W-1(т)Б^т. (7)

о

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

Хк+1 = Wk+^Хк + Б(/к, Хк )(/к+1 - ь)

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

Хк+1 = Wk+lWk Хк + (Б(/к, Хк) +

+Б(/к+1, Хк+1))

гк - гк-1 2 '

WJ■+1 = U(WJ■+1),

(6)

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

Хк+1 = Хк + Б(/к, Хк )(/к+1 - ч)

Используя формулу Симпсона, а также методы более высоких порядков, мы получим все многообразие подобных методов. Отметим только, что все эти методы решения интегрального уравнения (7) будут использовать приближенные значения фундаментальной матрицы ^ которые мы будем строить с учетом корректирующего правила (6). Изложим метод построения ортогональной матрицы полярного разложения матрицы W = Би, использующий аппарат матричных функций Ляпунова. Как было показано выше, в приложениях весьма важна задача построения полярного разложения невырожденной матрицы, т. е. представления ее в виде произведения симметричной положительно определенной матрицы на ортогональную. Такое представление всегда возможно и однозначно [2]. Пусть А — заданная невырожденная (п х п) -матрица. Требуется найти симметрическую положительно такие определенную матрицу Б и ортогональную матрицу и, что выполняется соотношение

4

Математические методы. Моделирование. Экспериментальные исследования

А = Ш. (8)

Введем в рассмотрение множество Мп квадратных матриц порядка п и множество Бп невырожденных матриц того же порядка. Множество Мп — линейное нормированное пространство с нормой, определяемой соотношением

Лсх||

|| С ||= тах

ХеЕп ||Х|| '

(9)

где С е Мп, Х е Еп. На множестве Мп можно задать динамическую систему, определяемую матричной системой дифференциальных урав-

нений

Х = Б(Х ^),

(10)

где Хе М ; для ¥ выполнены условия, обеспечивающие существование, единственность и продолжимость на интервале ж) решений; точка обозначает дифференцирование по параметру t.

Поставим следующую задачу: построить такую дифференциальную систему вида (10), что решение задачи Коши с начальным условием Х = Х0 при t = ^ = 0 сходится к значению матрицы и полярного разложения (8). Таким образом, мы сведем первоначальную задачу к численному интегрированию построенной системы дифференциальных уравнений. Рассмотрим матричное дифференциальное уравнение

Х = !(-Х + Х *_1). 2

(11)

Обозначим: Х ^, Х0) — решение системы (11), удовлетворяющее начальному условию Х = Х0 при t = 0 .

Теорема 1. Если матрица начальных условий Х0 невырожденная, то решение Х(^Х0) системы (11) также будет невырожденной матрицей при t > 0 , причем решение Х(1,Х0) неограниченно продолжимо при t > 0 [2].

Доказательство. Сначала покажем, что при t > 0 выполняется условие

Х(г,X,)еВп. (12)

Предположим противное, а именно: пусть существует момент времени ^, в который выполнено соотношение ёе1(Х(^,Х0)Х (t1,Х0))Ф 0. Тогда справедливо равенство

йе1(Х«ьХ0)Х*«ьХ0))Ф 0. (13)

Введем в рассмотрение матричную функцию

*

V = ХХ , заданную в Мп. Эта функция удовлетворяет на решениях системы (11) следующему уравнению:

V = _У + Е . (14)

Интегрируя, получаем [3]

V = (V - Е)ехрН) + Е. (15)

В силу того, что матрица % = Х0Х0 является симметрической и положительно определенной, из формулы (14) следует, что при t > 0 определитель матрицы Уне обращается в нуль, а это находится в противоречии с формулой (12), которая и доказывает справедливость условия (12). Необходимым и достаточным условием продолжимости решений системы (11) является рас-

ге .

ходимость интегралов у у аУу . Из (13) имеем

t

| у_ У у = _1п(уу _Ъу )|0 ,

0

где Ъу — символ Кронекера [4]. Расходимость интеграла следует из формулы (14), так как из нее вытекает у у ^ Ъу , и выражение в правой части последнего соотношения — неограниченное. Доказательство закончено.

Теорема 2. Решение Х(^Х0) системы (11) с начальным условием Х0 = А сходится при t к значению матрицы и полярного разложения (8), причем справедлива оценка

|| Х(/,А)_ и ||<|| АА* _ Е|| exp(_t). (16)

Доказательство. Поскольку решение Х^, А) удовлетворяет условиям теоремы 1, то оно пред-ставимо единственным образом в виде

Х ^ ,А) = F(t )U(t), (17)

где F(t), U(t) — матрицы полярного разложения. Отметим, что F(0) = F, и(0) = и. Подставляя выражение (16) в уравнение (11), можно получить следующее соотношение:

Б (/ ^) + F(t )1и (/ )и* (/ )F(t) =

= F(t)U(t)1и*(/)F(t) + F(t^(/). (18)

Из (14) следует, что матрица F(t) представляет собой ряд по целым отрицательным степеням экспоненциальной функции параметра t с коэффициентами, являющимися постоянными симметрическими матрицами:

F(t) = E + |(F2 - E)e- - -4(F2 - E)2 e~2t +...

Следовательно, имеет место тождество

F (t )F(t) = F(t )F (t). (19)

Рассматривая совместно выражения (17), (18) и условие ортогональности

U(t)U (t) + U(t)U (t) = 0, получаем, что имеет место соотношение [5]

U (t) = 0, или U(t) = const = U.

Поскольку || f (t)- E ||—> 0 в силу (14), то справедливо утверждение

||X(t,A)-U||— 0 при t —

Оценку (15) получаем следующим образом. По свойству нормы (19) можно выписать цепочку равенств и неравенств

|| Х(/, Л) - и || = || ) - Е)и(/) || < || ) - Е|| = = || (Б2 (г) - Е)(Б(/)+Е)-11| < || Б2 (/) - Е|| = = || Б2 - Е|| ехр(-) = || ЛЛ* - Е || ехр(г).

При оценке нормы || Б2(г) - Е|| мы использовали выражение (14). Доказательство закончено.

На основании доказанных теорем можно утверждать, что значение матрицы и можно получить, численно проинтегрировав систему (11) с начальным условием Х0 = Л при г = 0. Точность нахождения матрицы и будет зависеть как от длины интервала интегрирования, так и от точности самого метода интегрирования.

Работа выполнена при финансовой поддержке РФФИ (пр. № 10-08-000624).

СПИСОК ЛИТЕРАТУРЫ

1. Зубов, И.В. Анализ управляемых систем и равновесных движений [Текст] / И.В. Зубов, Н.В. Зубов, М.В. Стрекопытова / ВВМ/— СПб., 2012.— 322 с.

2. Зубова, А.Ф. Математические методы моделирования промышленных процессов и технологий [Текст] / А.Ф. Зубова / СПбГУ.— СПб., 2004.— 472 с.

3. Зубов, А.В. Математические методы исследования устойчивости и надежности технических систем [Текст] / А.В. Зубов, Н.В. Зубов, С.В. Зубов,

А.Ф. Зубова / ВВМ.— СПб., 2011.— 362 с.

4. Зубов, Н.В. Динамическая безопасность управляемых систем [Текст] / А.В. Зубов, Н.В. Зубов.— СПб.: Изд-во НИИ Химии СПбГУ, 2009.—172 с.

5. Зубов, С.В. Математические методы качественного анализа систем управления и устойчивость расчетных движений [Текст] / А.В. Зубов, С.В. Зубов.— СПб.: Изд-во АООТ «Мобильность-плюс», 2012.— 357 с.

REFERENCES

1. Zubov, I.V. The analysis of controlling systems and equally weight motions [Text] / I.V. Zubov, N.V. Zubov, M.V. Strecopitova / VVM.— SPb, 2012.— 322 p.

2. Zubova, A.F. The mathematics methods of modeling industry processes and technology [Text] / A.V. Zubova / SPbGU.— SPb., 2004.— 472 p.

3. Zubov, A.V. The mathematics methods of investigation of stability and safety of technical systems [Text] /

A.V. Zubov, N.V. Zubov, S.V. Zubov, A.F Zubova / VVM.— SPb., 2011.— 362 p.

4. Zubov, N.V. The dynamics safety of controlling systems [Text] / A.V. Zubov, N.V. Zubov.— SPb.: Publisher NII of Chemistry SPbGU, 2009.— 172 p.

5. Zubov, S.V. The mathematics methods of qualitative analysis of systems control and stability of calculating motions [Text] / A.V. Zubov, S.V. Zubov— SPb.: AOOT «Mobility-plus», 2012.— 357 p.

СВЕДЕНИЯ ОБ АВТОРАХ /AUTHORS

ЗУБОВ Афанасий Владимирович — доктор физико-математических наук профессор факультета прикладной математики — процессов управления Санкт-Петербургского государственного университета; 198504, Старый Петергоф, Университетский пр., д. 35, г. Санкт-Петербург, Россия; e-mail: [email protected]

ZUBOV Afanasiy V. — St.-Petersburg State University; 198504, Stariy Petergof, Universitetski'y pr., 35, St.-Petersburg, Russia; e-mail: [email protected]

Математические методы. Моделирование. Экспериментальные исследования -►

КАЛЯДА Леонид Герасимович — аспирант факультета прикладной математики — процессов управления Санкт-Петербургского государственного университета; 198504, Старый Петергоф, Университетский пр., д. 35, г. Санкт-Петербург, Россия; e-mail: [email protected]

KALYADA Leonid G. — St.-Petersburg State University; 198504, Stariy Petergof, Universitetski'y pr., 35, St.-Petersburg, Russia; e-mail: [email protected]

НЕЧАЕВ Алексей Иванович — аспирант факультета прикладной математики — процессов управления Санкт-Петербургского государственного университета; 198504, Старый Петергоф, Университетский пр., д. 35, г. Санкт-Петербург, Россия; e-mail: [email protected]

NECHAEV Alexei I. — St.-Petersburg State University; 198504, Stariy Petergof, Universitetski'y pr., 35, St.-Petersburg, Russia; e-mail: [email protected]

УЖЕГОВ Николай Сергеевич — аспирант факультета прикладной математики — процессов управления Санкт-Петербургского государственного университета; 198504, Старый Петергоф, Университетский пр., д. 35, г. Санкт-Петербург, Россия; e-mail: [email protected]

UGEGOV Nicolay S. — St.-Petersburg State University; 198504, Stariy Petergof, Universitetski'y pr., 35, St.-Petersburg, Russia; e-mail: [email protected]

© Санкт-Петербургский государственный политехнический университет, 2013