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

Математика. Физика

УДК 519.622.2

О ПОСТРОЕНИИ ПЕРИОДИЧЕСКИХ РЕШЕНИЙ ОДНОГО КЛАССА НЕАВТОНОМНЫХ СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ В РАСПРЕДЕЛЕННОЙ КОМПЬЮТЕРНОЙ СРЕДЕ

А.Н. Пчелинцев1, А.Ю. Поветьев2, В.Е. Подольский2

Кафедры: «Прикладная математика и информатика» (1),

«Системы автоматизированного проектирования» (2), ГОУ ВПО «ТГТУ»; [email protected]

Представлена членом редколлегии профессором В. И. Коноваловым

Ключевые слова и фразы: периодическое решение; ряд Фурье; символьные вычисления; система обыкновенных дифференциальных уравнений; схема последовательных приближений.

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

Введение

Довольно часто (и на практике в том числе) возникает задача построения периодических решений нормальной системы обыкновенных дифференциальных уравнений вида

X = f (^ x), (1)

где x(t) - векторная функция действительного переменного Ґ, f - векторная функция, равная

f (^ x) = ф( x) + h(t), (2)

где вектор ф(x) - многомерный многочлен, а функция h(t) представляет собой векторнозначный тригонометрический полином (Г-периодическая векторная функция).

Многие из теорем существования периодических решений у системы (1)

(см., например, [1]) используют тот фундаментальный факт, что такие решения

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

Пусть известно, что система (1) имеет единственное Т-периодическое реше-

*

ние x (t). Примерами систем, имеющих единственное периодическое решение, являются системы с конвергенцией [2, 3]. Рассмотрим один класс таких систем,

*

для которого предложен способ построения приближения к решению x (:); приведем условия, накладываемые на функцию f позволяющие показать сходимость схемы последовательных приближений на периоде. При этом введем вспомогательную систему для построения приближений в символьном виде к некоторой периодической функции, зависящей от начальных условий для системы (1). Варь*

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

1. Условия, накладываемые на систему обыкновенных дифференциальных уравнений вида (1)

Будем рассматривать класс систем (1), подчиняющийся следующим условиям.

1. Система (1) является системой с конвергенцией, то есть она имеет единст-

*

венное Т-периодическое решение x (О, асимптотически устойчивое в целом.

*

2. Замкнутый шар Бг радиуса г, который содержит значения функции х (:),

лежит внутри замкнутого шара радиуса Я. Для векторов р и q из £я имеет

место неравенство Липшица

I ф( р) -ф^)1 < 1I р - q I, (3)

где положительное число I удовлетворяет условию

I < —. (4)

2Т

3. Найдется такое положительное число М, что для всех векторов х из £я

| х) | < М, г + 2ТМ < Я. (5)

2. Переход к вспомогательной системе уравнений

Для простоты обозначений положим, что начальный момент времени равен нулю. Перепишем (1) в интегральной форме

г

х^) = С + | f (т, х(т)) dт,

о

*

где вектор С определяет начальные условия. Поскольку х (:) - Т-периодическая функция, то

х* (0) = х* (Т) = С *.

Тогда

Т

I f (т, (* (т))dт = 0. (6)

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

Поскольку функция х (t) удовлетворяет системе (1), то она удовлетворяет и системе

1

у = f (t, У) - — J f (Т У(т)) dT. T 0

(7)

ds.

(8)

Перейдем от равенства (7) к интегральному соотношению

t Г 1 т

y(t) = С + | /(5, y(s)) - — j /(т, y(x))dx

0 L 0 _

Заметим, что переход к вспомогательной системе (8) необходим, поскольку вложенный интеграл (среднее интегральное значение на периоде) позволяет при вычислениях избежать появления в символьных выражениях степеней t. Это уменьшает объем памяти, отводимой под них, а также формирует символьное представление искомой функции как приближение к разложению в ряд Фурье периодического решения системы (1).

3. Схема последовательных приближений

*

Чтобы получить приближение к решению х (t) на периоде, сначала построим приближение к функции, являющейся решением системы (8). Затем, варьируя в ней вектором С, найдем требуемое приближение. Для этого покажем, что мы находимся в условиях применимости метода последовательных приближений.

Пусть fi - пространство непрерывных Г-периодических векторных функций y(t) со значениями из шара SR . Расстояние между функциями p,q efi определим как

р(p, q)= max | p(t) - q(t)|.

te[0,T ]

Таким образом, пространство fi является метрическим. В fi рассмотрим оператор

t 1 T

Фу (t) = C + J f (s, y(s)) - T J f (t, У(т)) dx

t

ds = C + J u(s)ds. 0

(9)

Покажем, что функция

ё (0 = Фу(0 принадлежит пространству О , когда у еП.

*

В силу (9) и | С | < г (мы будем искать вектор С в шаре Бг, поскольку в нем

*

содержатся значения функции х )) имеем

|g (t)|<|C| +

J u(s)ds

< Г + 2T max |f (t, y(t ))|.

te[0,T ], yeQ

Учитывая, что y(t) принимает значения из £я , и неравенства (5), получим

|ё (t)| < Я

для всех t е [0,Т].

Теперь покажем Т-периодичность функции ё (г)

ё (t + Т) - ё (t ) = 0

для любых действительных значений t. С учетом (9) рассмотрим разность

t+T t

g (t + T ) - g (t )= | u(s)ds -J u(s)ds = | u(s)ds + J u(s)ds. 0 0

t+T

= Ju(

0

Поскольку функция f Т-периодична по г и у (г) также Т-периодична, то функция и^) будет Т-периодической. С учетом этого, делая замену s = а + Т в первом интеграле, получим

г оо Т

ё (г + Т) - ё (г )= | u(a)da + | и^^ = | и^^ = | u(s)ds.

-т

Положим

Тогда

B = J f (т, y(x))dx.

g (t + T ) - g (t) = J

f (s, y(s)) - jB

ds = B-------BT = 0.

T

Таким образом, если у еО, то ё еП .

Покажем, что отображение Фу сжатое. Оценим

р(Фр, Фд)= max

te[0,T ]

ds

I f (s, p(s)) - f (s, q(s)) - Т IУ(т, р(т)) - f (т, q(т))}dт

о Т о

В силу неравенства (3) и соотношения (2), получаем

р(Фр, Фq) < 21Тр(р, q).

Из (4) следует, что

21Т = I <1.

Таким образом, отображение является сжатием. Теперь покажем полноту пространства О.

Следуя формуле (9), образуем бесконечную последовательность

Ут = ФУт-1.

элементов уо, у1, У2, ... из О. Поскольку оператор Ф сжимающий, запишем последовательность неравенств:

р( У2, уО < ^р( Уl, уо);

рО^У2) ^LP(У2,у{) ^L р(уъУ0);

Р(Ут+1, Ут ) ^ LmР(Уь У0).

По неравенству треугольника имеем

р(ут+], ут ) < р(ут+], ут+]-1) + р(ут+]-1, ут+]-2 ) + +

+ Р(ут+1, ут ) < (V -1 + Ь -2 + • • • + Ь +1)Ьтр(у1, уо) <

(1 - Ь) т , ч Ьт ,

< ——— Ь P(У1,уо) <-—гР(уъуо)<^

1- Ь 1- Ь

когда т > Мв, где Мв достаточно велико. Седовательно, указанная последовательность является последовательностью Коши по определению. Поскольку

| ут (г) < Я

то по теореме Больцано-Вейерштрасса из данной поледовательности можно выделить сходящуюся подпоследовательность. Таким образом, по следствию 2 [4, с. Ю1] последовательность Коши также будет сходиться в пространстве О. Поэтому это пространство является полным.

Тогда, согласно методу последовательных приближений, схема

Т

>m+1(t) - C+J

f (А Ут (s)) - T J f (T Ут (т))т

T 0

ds. (10)

0 L

сходится к Г-периодическому решению системы (8).

4. Нахождение периодического решения

Поскольку правая часть системы (1) по х является многомерным многочленом, а по t представляет собой тригонометрический полином, начальную функцию y0 (t) целесообразно выбирать равной вектору С или как

y0 (t ) = С cos(rat),

где ю = 2п/Т, то есть y0 e fi.

Заметим, что произведение двух (а, следовательно, и любого числа) тригонометрических функций (cos или sin) кратных аргументов, а также степени таких функций (см., например, [5, с. 39]), можно представить в виде суммы постоянного вектора и тригонометрического полинома, и при интегрировании тригонометрического полинома получается такая же сумма. Тогда из схемы (10) каждую итерацию можно вычислять символьно. При этом преобразование тригонометрических функций в символьном виде, а также их символьное интегрирование, распараллеливается. Идея параллелизма заложена в формуле (10): вычисления для каждой компоненты вектора ym+i(t) можно производить независимо друг от друга, а полученные символьные выражения хранить в сетевой базе данных, доступной считающим процессам в распределенной компьютерной среде. Здесь база данных представляет собой структурированную общую память для параллельных процессов.

Следуя формуле (10), мы строим функцию ym+i(t) до такого значения m, когда

max р(ym+1, ym ) < e p,

CeSr

где e p - точность для схемы (10).

Рассмотрим интеграл

T

I = j f (т, y(T))dx.

0

При подстановке в него решения y(t) системы (8) получим векторную функцию от С. Тогда переобозначим

5(C ) = I.

Для того чтобы схема (10) сходилась к периодическому решению системы (1),

*

в силу (7) нужно подобрать вектор C так (C = C ), чтобы имело место равенство (6), то есть

5(с *)= 0. (11)

Поскольку система (1) имеет единственное периодическое решение, то система алгебраических уравнений (11) в области Sr будет также иметь единственное решение.

Система (11) равносильна тому, что

(5(с'),5(с' ty=°.

Откуда получаем задачу оптимизации для нахождения приближенного значения

*

вектора C

0(C) = (5(C), 5(C)) ^ min, C e Sr. (12)

Переход от системы (11) алгебраических уравнений к задаче (12) связан с тем, что непосредственно при расчетах функцию y(t) мы определяем приближенно.

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

Если порядок велик, то в этом случае целесообразно записать теорему Кару-ша-Джона (см., например, [6, с. 240]), что сведет решение задачи (12) к решению системы алгебраических уравнений. Если возможно здесь применение метода Ньютона, то вычислительный процесс также можно распараллелить.

Метод Ньютона включает в себя нахождение матрицы, обратной матрице Гессе. Это связано с большими вычислительными затратами. Однако как процесс нахождения самой матрицы Гессе, так и ее обращение параллелится. Заметим, что компоненты (значения частных производных второго порядка и смешанных производных) матрицы Гессе целевой функции могут быть определены независимо друг от друга. В этом случае данная задача распределяется по процессам, каждый из которых вычисляет соответствующую производную. Используя теорему Гамиль-тона-Кэли для матрицы Гессе, можно получить выражение для ее обратной матрицы [7], вычисляя которую потребуется определить степени исходной матрицы. Здесь процесс умножения матрицы на матрицу также распараллеливается [8].

5. Пример

В качестве примера использования описанного способа рассмотрим нелинейное уравнение колебаний (в данном пункте мы переобозначим некоторые функции, как принято в известной литературе)

X + q( x) X + g (x) = p(t), (13)

где q(x) = 0,01(1 + 2x2), g(x) = 0,03x, p(t) = 0,01cos(24t), период T правой части уравнения (13) равен п/12. По теореме 8.1 [2, с. 108] для уравнения (13) выполняется свойство конвергенции.

От уравнения (13) перейдем к эквивалентной системе второго порядка [2, с. 107, 110], правая часть которой имеет вид (2)

x = y - Q(x) + P(t), y = -g (x), (14)

где

x ( 2 ^ t 1 Q( x) = J q(Qd^ = 0,01^ x + - x3 J, P(t) = J p(x)dx = -^^sin^t).

Найдем радиус r шара Sr, в которую погружаются все решения системы (14). Для этого сначала покажем выполнение условий 1-4 [2, с. 109].

1. Функции q, g и p непрерывны, g удовлетворяет условию Липшица с константой 0,03.

2. Значения чисел а, а ив выберем равными 1/12, 73/7200 и 1/400 соответственно, что q(x) > а при | x | > a, g(x) > в при x > a, g(x) < -в при x < -a.

3. Функция P(t) ограничена при всех t, то есть | P(t) |< E, E = 1/1200.

4. Существует такое число у = 1/259200, что Q(x) - E >у при x > a ,

Q(x) + E < -у при x < -a ; G(x) > 0 при | x |> a, где

x

G( x) = J g (§)d| = 0,015x2.

0

Тогда любая траектория системы (14) навсегда погружается в прямоугольник, задаваемый неравенствами

| x | < a, | y | < Х0 + X = b,

где ^0 таково, что | Q(x) | +E < X0 при | x | < a. Выберим его равным

^0 = max | Q(x) | +E.

|x|<a

Поскольку Q(x) - всюду возрастающая функция, то X0 = у . Значение постоянной X выбирается как [2, с. 113]

[ 8 „

X1 — |A0a—+ А,

V Y '

где А - произвольное положительное число и

8 = max | g (x) |

| x|< a

[2, с. 110]. Положим А = >/3/1200-у . Тогда b = V3/600 .

Таким образом, радиус r можно принять равным

I 2 + b2 V7

r = V a + b =------.

30

Теперь установим выполнение условий 2 и 3 из п. 1. Перепишем систему (14) в векторной форме

Z = w(t, z),

где

x(t ) " y — Q( x) "1/2400"

, w(t, z) = ф(z) + h(t), ф^) = , h(t ) = sin(24t)

_ y(t )_ _ — g ( x) _ 0

z(t ) =

Будем рассматривать шар , внутри которого находится шар 8Г, радиуса Я = 0,5. Тогда из (14)

г + 2Т тах | м>Ц, г($)) | <

tе[0,T], геО

< r + 2T max

|x|<R, |y|<R

I y I +0,01| |r|+ f|x|3| + * + 0.031 XI

< R.

Локальную липшицевость для векторной функции ф(г) установим следующим образом. Представим

ф(2) - ф() = А(XI, Х2)(21 - г2),

где матрица A имеет вид

A( xh x2) =

Отсюда получаем

— 0,01| 1 + "з {xi + X1X2 + X2 } I 1

- 0,03

0

2T max Щ( xi, X2 )| < 1.

|x^<R. |x2 |<R

6. Вычислительный эксперимент

Для поиска приближения к периодическому решению системы (14) был разработан комплекс программ для кластера Тамбовского государственного технического университета. Данный комплекс составляют две программы:

1) построение символьных выражений по формуле (10);

2) поиск экстремума функции 9n (С), являющегося приближением к начальному значению периодического решения системы (14), для итерации с номером n.

Первая программа на каждой итерации вызывает математический пакет Maxima символьных вычислений для разложения степеней тригонометрических функций в тригонометрический полином (функция пакета trigreduce) и взятия от полученного выражения неопределенного интеграла (функция integrate). Заметим, что после интегрирования из символьного выражения исключаются члены, содержащие линейную функцию от t вида kt (в силу наличия в формуле (10) среднеинтегрального значения).

Поскольку система (14) нелинейна, от итерации к итерации вычисляемые символьные выражения резко увеличиваются в размерах. Поэтому структура таблицы Fourier базы данных на языке SQL выглядит следующим образом:

CREATE TABLE Fourier (

ind TEXT, arg TEXT, koef LONGTEXT

);

Здесь поле ind хранит в себе номер итерации и номер координаты, поле arg хранит тригонометрическую функцию с аргументом или 1 (для свободного члена) и поле

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

В качестве начальной функции была взята вектор-константа C. На первой и

второй итерациях обнаружено, что значения найденного вектора C* для каждой из них одинаковы. Функция 02(C) имеет вид (данные, полученные с использованием пакета Maxima)

((

hC ) =------

2 144

c2

c1

c1

477757440000000000 5503765708800000000000

169075682574335946915840173c

507227047723007840747520199c

25361352386150400000000000000 1460813897442263040000000000000000

68155733137438134529454363209482230189c

6815573319906622439424000000000000000000

42

68155733164828395106496786609848317479 392577023226621452510822400000000000000000000

+ 1^ + -

1

100 1920000

2

где

C =

Цикл, соответствующий найденному начальному значению С

2

+

+

х

Откуда значение вектора С , определяющего искомое приближение к решению

*

х (/),равно

^-1/57600"

0

С * =

Полученный вектор C был проверен с помощью комплекса программ, описанного в работе [9], на возврат траектории системы (14) в окрестность начального значения через период (рисунок).

Работа выполнена при поддержке РФФИ (проекты 10-07-00136 и 11-0700098).

Список литературы

1. Красносельский, М. А. Оператор сдвига по траекториям дифференциальных уравнений / М. А. Красносельский. - М. : Наука, 1966. - 332 с.

2. Плисс, В.А. Нелокальные проблемы теории колебаний / В.А. Плисс. -М. ; Л. : Наука, 1964. - 367 с.

3. Демидович, Б.П. Лекции по математической теории устойчивости / Б.П. Демидович. - М. : Наука, 1967. - 472 с.

4. Шварц, Л. Анализ. В 2 т. Т. 1 / Л. Шварц ; пер. Б.П. Пугачева ; под ред. С.Г. Крейна. - М. : Мир, 1972. - 824 с.

5. Градштейн, И. С. Таблицы интегралов, сумм, рядов и произведений / И.С. Градштейн, И.М. Рыжик. - М. : Физматгиз, 1963. - 1100 с.

6. Поляк Б.Т. Введение в оптимизацию / Б.Т. Поляк. - М. : Наука, 1983. - 384 с.

7. Пчелинцев, А.Н. Определение параметров функции управления при решении задачи энергосберегающего управления в распределенной компьютерной среде / А.Н. Пчелинцев, В.А. Погонин // Системы упр. и информ. технологии. -2010. - № 1 (39). - С. 46-49.

8. Шпаковский, Г.И. Программирование для многопроцессорных систем в стандарте MPI / Г.И. Шпаковский, Н.В. Серикова. - Минск : БГУ, 2002. - 323 с.

9. Пчелинцев, А.Н. Численные методы построения обобщенно-периодических решений дифференциальных уравнений при моделировании динамических процессов : дис. ... канд. физ.-мат. наук : 05.13.18 : защищена 12.11.2009 : утв. 12.02.2010 / Пчелинцев Александр Николаевич. - Воронеж, 2009. - 114 с.

About the Construction of Periodical Solutions of One Class of Non-Autonomous Systems of Differential Equations in a Distributed Computing Environment

A.N. Pchelintsev1, A.Yu. Povetyev2, V.E. Podolskiy2

Departments: «AppliedMathematics and Computer Science» (1), «Computer-AidedDesign» (2), TSTU; [email protected]

Key words and phrases: Fourier series; periodic solution; the scheme of successive approximations; symbolic computation; the system of ordinary differential equations.

Abstract: The paper proposes the technique for the construction of approximations to the periodical solutions of one class of non-autonomous systems of common differential equations. It is based on the scheme of successive approximations with the use of symbolic calculations in a distributed computing environment to obtain the solution in an analytic form. The fitting of the scheme of successive approximations in the period is shown and also the example of the second grade system with the described calculation scheme is considered. The results of computational experiment are presented.

Über Konstruktion der periodischen Lösungen der einen Klasse der nichtautonomen Systemen der Differentialgleichungen im verteilten Computermedium

Zusammenfassung: Es ist das Verfahren des Baues der Approximationen zu den periodischen Lösungen der einen Klasse der nichtautonomen Systemen der gewöhnlichen Differentialgleichungen vorgeschlagen. Zugrunde ist das Schema der aufeinanderfolgenden Approximationen mit der Benutzung der Symboberechnungen im verteilten Computermedium für das Erhalten der Lösung in der analytischen Form gelegt. Es sind die Resultate des Berechnungsexperimentes angeführt.

Sur la construction des solutions périodiques d’une classe des systèmes non autonomes des équations différentielles dans un milieu informatique distributive

Résumé: Est proposé un moyen de la construction des approximations vers les solutions périodiques d’une classe des systèmes non autonomes des équations différentielles ordinaires. A la base est mis le schéma des approximations consécutives avec l’emploi des calculs symboliques dans un milieu informatique distributif pour l’obtention des solutions dans une forme analytique. Est montrée la convergence du schéma des approximations consécutives sur une période ainsi qu’est examiné un exemple du système du deuxième ordre pour lequel est appliqué le schéma indiqué des calculs. Sont cités les résultats de l’expériment de calcul.

Авторы: Пчелинцев Александр Николаевич - кандидат физико-математических наук, доцент кафедры «Прикладная математика и информатика»; Поветьев Алексей Юрьевич - аспирант кафедры «Системы автоматизированного проектирования»; Подольский Владимир Ефимович - доктор технических наук, профессор, проректор по информатизации, заведующий кафедрой «Системы автоматизированного проектирования», ГОУ ВПО «ТГТУ».

Рецензент: Туголуков Евгений Николаевич - доктор технических наук, профессор кафедры «Техника и технологии производства нанопродуктов», ГОУ ВПО «ТГТУ».