Методы анализа устойчивости экономических дифференциальных моделей

5. Москалева, Р. Н. Реализация принципа преемственности в обучении учащихся начальной и основной ступеней школы с углубленным изучением математики Автореф. дис. ... канд. пед. наук - Магнитогорск, 2007. - 24 с.

6. Стойлова, Л. П. Способы решения стохастических задач // Начальная школа. 1994. - №1. - С.84-86.

7. Селютин, В.Д. Научные основы методической готовности учителя математики к обучению школьников стохастике: Дис. ... д-ра пед. наук. - М., 2002. -344 с.

8. Цулина, И. В. Методические особенности курса «Элементы теории вероятности» для учащихся старших классов, Автореф. дис. ... канд. пед. наук, - М., 2010. - 20 с.

УДК 517.91: 518.1 ББК 22.193

В.Н. Кузнецов., О.Н. Папко

МЕТОДЫ АНАЛИЗА УСТОЙЧИВОСТИ ЭКОНОМИЧЕСКИХ ДИФФЕРЕНЦИАЛЬНЫХ МОДЕЛЕЙ

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

Ключевые слова: Дифференциальные модели экономики, исследование устойчивости, разностные схемы.

V.N. Kuznetsov., O.N. Papko

METHODS OF ANALYSIS OF ECONOMIC STABILITY DIFFERENTIAL MODELS

Absrtact. The article provides an overview of methods for analyzing the stability of systems of differential equations. Given the parameters are performed stability analysis of the Lorenz system, describing the dynamics of small urban systems are part of the metropolis.

Key words: Differential model of the economy, the study of stability, finite difference schemes.

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

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

^ = аг,

где А =

dt

- квадратная матрица постоянных коэффициентов.

Пусть Р(X ) = а0Хп + а1Хп +...+ап - некоторый многочлен, причём а0, а1,...,ап - действительные числа и а 0 > 0. Образуем матрицу п х п :

а1 а0 0 ... 0

M =

0 0 0 ... ап

По главной диагонали этой матрицы откладываются коэффициенты

ке от этих элементов расположены коэффициенты с убывающими номерами, влево - с возрастающими. При этом полагается а г = 0 , если г < 0 или г > п . Такая матрица М называется матрицей Гурвица. Главные диагональные миноры этой матрицы будут иметь вид

г. Вправо по стро-

a

a

1 n

а

а

n1

nn

a3 a2 a1 ... 0

a5 a4 a3 ... 0

Aj = al , A2 =

4 "0

a о a ^

A«=a»A B-i.

«3 "2

Критерий Гурвица. Для того, чтобы стандартный полином P (X ) был полиномом Гурвица, необходимо и достаточно, чтобы все главные диагональные миноры его матрицы Гурвица были положительными

Ak > 0, k = 1,...,n .

Выполнение критерия Гурвица является необходимым и достаточным условием устойчивости решения системы линейных ОДУ с постоянной матрицей коэффициентов [1].

Устойчивость систем линейных ОДУ с почти постоянной матрицей. Если система (1) устойчива при t ^ да, то система

= [A + B

dt

где B (t) - матрица переменных коэффициентов, для которой выполнены следующие ограничения

1) каждый элемент матрицы B (t) есть непрерывная функция на промежутке [t0,<»),

да

2) ЛB(t)|| dt<да

t0

также устойчива при t ^ да .

Если система (1) асимптотически устойчива при t ^да, то возмущённая линейная система (2) также асимптотически устойчива, если для матрицы B (t) выполняется условие 1) и B (t) ^ 0 при

t ^ да .

Характеристический показатель Ляпунова. Достаточное условие асимптотической устойчивости. Число, определяемое формулой

Х[ t ]=iimiin| f (t )|,

t^да t

называется характеристическим показателем Ляпунова. Если матрица линейной системы ограничена

то каждое нетривиальное решение Y = Y (t) имеет конечный характеристический показатель. Множество всех характеристических показателей решений системы называется её спектром. Для асимптотической устойчивости однородной линейной системы достаточно [1], чтобы наибольший её характеристический показатель был отрицательным

а = max ak < 0 .

k

Приводимые системы. Теорема Н.П. Еругина. Пусть снова рассматривается система линейных ОДУ.

Матрица L (t) непрерывная и непрерывно дифференцируемая на промежутке [t0,<x>) называется матрицей Ляпунова, если

1) L(t) и L'(t) ограничены на промежутке [t0,<»).

2) | det L (t )| > m , m = const, m > 0. Линейная система

^ = A

dt

называется приводимой, если существует такая матрица Ляпунова L(t) , что подстановка

X = L (t )Y

приводит систему к системе уравнений

— = BX, dt

где B - постоянная матрица.

Теорема Еругина (необходимое и достаточное условие устойчивости). Система линейных ОДУ (3) приводима тогда и только тогда, когда некоторая её фундаментальная матрица Y (t) может быть представлена в виде матрицы Ляпунова L (t), умноженной на экспоненциал произведения постоянной матрицы B на независимую переменную t

Y(г) = L(г)eBt.

Для приводимой системы справедливы следующие условия устойчивости.

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

2) Приводимая система линейных ОДУ асимптотически устойчива тогда и только тогда, когда все её характеристические показатели отрицательны [1].

Системы линейных ОДУ с периодическими коэффициентами. Если для системы (3) выполняется условие A( t + T ) = A(t ), то система имеет по крайней мере одно решение Y ( t ), не равное тождественно нулю, такое, что

Y (t + T ) = X Y ( t ),

для всех t, X = const, X* 0 [1].

Для линейной системы (3) с T - периодической матрицей нормирования при t = 0 фундаментальная матрица решений имеет вид

Y ( t ) = Ф ( t ) e Л t ,

где Ф ( t ) - непрерывная и непрерывно дифференцируемая T - периодическая неособенная матрица, Ф (0) = E, Л - постоянная матрица. Матрица Y (T) называется матрицей монодромии.

Собственные значения X j матрицы Л называются характеристическими показателями системы (1). _

Собственные значения pj, ( j = 1,..., n ) матрицы C = Y ( T ) называются мультипликаторами. Однородная линейная периодическая система с непрерывной матрицей устойчива тогда и только тогда, когда все её мультипликаторы pj расположены внутри замкнутого единичного круга | p | < 1

, причём мультипликаторы, лежащие на окружности | p | = 1, имеют простые элементарные делители, если их рассматривать как собственные значения соответствующей матрицы монодромии. Для асимптотической устойчивости периодической системы необходимо и достаточно, чтобы все её мультипликаторы находились внутри единичного круга | p | < 1.

Второй метод Ляпунова. Второй (качественный, прямой) метод Ляпунова в основном связан с введением вспомогательных оценивающих функций, получивших его имя, и явился, по словам историков науки, результатом синтеза идей качественной теории дифференциальных уравнений и исключительной интуиции А. Пуанкаре с аналитическим гением А.М. Ляпунова [1]. Метод функций Ляпунова стал основным методом классической и современной теории устойчивости и качественной теории дифференциальных уравнений исключительно эффективным для разнообразных нелинейных систем, особенно не являющихся квазилинейными [1]. Пусть дана система уравнений

dy,

= fj

dt

Предполагается, что функции f (t,y1,...,yn ) непрерывны в некоторой открытой области Z , которая может совпадать со всем пространством. Кроме того функции f (t,y1,...,yn ) удовлетворяют в любой замкнутой области G , лежащей в Z , условиям Липшица. В этих условиях существует единственное решение системы (4), удовлетворяющее заданным начальным условиям [1]. Функция V(t, Y) ,определённая в фазовом пространстве переменных y1, ...,yn , непрерывная в некоторой области Z , включающей в себя начало координат, называется определённо положительной в области Z, если всюду в области Z, кроме точки O(0,...,0) имеет место неравенство V (t ,Y ) > 0.

Если же выполняется неравенство V (t ,Y ) < 0, то функция V( t, Y ) называется определённо отрицательной.

Если в области Z имеет место всюду неравенство V(t,Y)>0 или V(t,Y)<0, то функция V(t,Y) называется знакопостоянной.

В первом случае - знакоположительной, во втором - знакоотрицательной.

Теорема Ляпунова об устойчивости. Если для системы (4) существует в области Z знакоопреде-лённая функция V(t, Y), производная которой по времени V(t, Y), взятая в силу системы (4), является знакопостоянной функцией знака, противоположного знака функции V( t, Y ), то тривиальное решение системы (4) устойчиво в смысле Ляпунова.

Теорема Ляпунова об асимптотической устойчивости. Если для системы (4) существует знакооп-ределённая функция V( t, Y ), полная производная которой по времени, найденная в силу системы

(4), будет также знакоопределённой, знака противоположного с V(t, Y), то тривиальное решение системы (4) будет асимптотически устойчиво.

Теоремы о неустойчивости. Первые две теоремы принадлежат А.М. Ляпунову, третья теорема доказана Красовским [1].

Теорема 1. Если существует функция V(t, Y), имеющая знакоопределённую производную по времени, и такая, что в любой окрестности точки O функция V(t, Y) не является знакопостоянной, знака противоположного с V (t, Y), то нулевое решение системы (4) неустойчиво. Теорема 2. Если существует функция V (t, Y) такая, что её производная по времени имеет вид

dV = XV + W , dt

где X - положительная постоянная, а W или тождественно обращается в нуль или является знакопостоянной, и если в последнем случае функция V (t, Y) не является в любой окрестности точки O знакопостоянной, знака противоположного с W , то нулевое решение системы (4) неустойчиво. Теорема 3. Если существует функция V(t, Y), не являющаяся знакоотрицательной в произвольной окрестности точки O , и такая, что

dV dV

— > 0 вне M, — > 0 на M, dt dt

где M - множество, не содержащее целых траекторий ( кроме точки O ), то нулевое решение неустойчиво.

Способы построения функций Ляпунова. Первые способы построения функций V(t, Y) были предложены А.М. Ляпуновым [1]. Классические способы построения функций Ляпунова подразделяются на следующие.

a) Использование в качестве функций Ляпунова первых интегралов или их связки.

b) Использование структуры производной функции Ляпунова и решение матричного алгебраического ATB + BA = -C (для системы дифференциальных уравнений с постоянной матрицей A ли-

dB

нейной части) или матричного дифференциального--+ AT(t)B + BA(t) = -C ( для системы диф-

dt

ференциальных уравнений с матрицей A линейной части, зависящей от t ) уравнений Ляпунова. Здесь C - заданная постоянная симметрическая положительно определённая матрица, а B - искомая симметрическая матрица, определяющая функцию Ляпунова в виде квадратичной формы

V(Y) = YTBY [1].

c) Способы построения функций Ляпунова, основанные на задании их структуры в виде: интеграла энергии для механических и электромеханических систем [1]; псевдоквадратичных форм V(t,Y) = YTB(t, Y)Y;

линейных форм модулей V (t ,Y) = £ a ^x^ .

Для проверки условий положительной определённости функций Ляпунова и отрицательной определённости её производной использовались условия Рауса и Рауса - Гурвица. Ниже на основе метода анализа устойчивости, представленного в [2], выполняется анализ устойчивости решения модели Лоренца

dy1

— = стy2 -СТУ1, dt

dy2

— = -У2 + ГУ1 - У1У3 j (5) dt

dy3 ,

— = -ЬУ3 + УьУ2> dt

8

для заданного набора параметров - ст = 40, b = —, r = 80.

Критерий сводится к проверке неравенства (6). В случае его выполнения предполагается, что решение системы (4) устойчиво [2].

Ук

yk 0 - yk 0

< с , с = const, Vte [t0, да) , k = 1,...,n .

Программная модель, реализующая критерий (6) имеет вид:

program ust_ejler_lorenz; {$APPTYPE CONSOLE}

uses SysUtils; const

h=0.00001; n=3; t0=0; TT=10000; sigma=4 0; r=80; b=8/3;

eps1=0.0001; eps2=0.0001; eps3=0.0001; var

t, norma, y1, y2, y3, yv1, yv2, yv3, y11, y22, yv11, yv22 : extended; k, k0 : longint;

delta : array[1..n] of extended; yv : array[1..n] of extended; y01, y02, y03, yv01, yv02, yv03 : extended; function f1(t, y1, y2, y3 : extended) : extended;

begin f1:=sigma*(y2-y1); end; function f2(t, y1, y2, y3 : extended) : extended;

begin f2:=-y2+(r*y1)-(y1*y3); end; function f3(t, y1, y2, y3 : extended) : extended;

begin f3:=(-b*y3)+(y1*y2); end; begin

y01:=0; y02:=0; y03:=0;

yv01:=y01+eps1; yv02:=y02+eps2; yv03:=y03+eps3; t:=t0; k:=0; k0:=0; y1:=y01; y2:=y02; y3:=y03; yv1:=yv01; yv2:=yv02; yv3:=yv03;

delta[1]:=yv01-y01; delta[2]:=yv02-y02; delta[3]:=yv03-y03; repeat

y11:=y1; y22:=y2; yv11:=yv1; yv22:=yv2; y1:=y1+h*f1(t, y1, y2, y3); y2:=y2+h*f2(t, y11, y2, y3); y3:=y3+h*f3(t, y11, y22, y3); yv1:=yv1+h*f1(t, yv1, yv2, yv3); yv2:=yv2+h*f2(t, yv11, yv2, yv3); yv3:=yv3+h*f3(t, yv11, yv22, yv3); yv[1]:=yv1-y1; yv[2]:=yv2-y2; yv[3]:=yv3-y3; if ( abs( yv[1]/delta[1] )>=abs(yv[2]/delta[2] ( abs(yv[1]/delta[1])>=abs(yv[3]/delta[3] then norma:=abs(yv[1]/delta[1]); if ( abs(yv[2]/delta[2])>=abs(yv[1]/delta[1]) ( abs(yv[2]/delta[2])>=abs(yv[3]/delta[3]) then norma:=abs(yv[2]/delta[2]); if ( abs(yv[3]/delta[3])>=abs(yv[1]/delta[1] ) ( abs(yv[3]/delta[3])>=abs(yv[2]/delta[2]) then norma:=abs(yv[3]/delta[3]); k:=k+1; k0:=k0+1; t:=t+h;

if k>=50000000 then

begin writeln('t=',t:4:0,' ','norma=',norma) until t>=TT; readln; readln; end.

Результаты компьютерного моделирования системы (5) представлены в таблице 1.

Таблица 1

Значения нормы по критерию

t norma t norma

500 5.81896556134513E+0005 5500 9.65082710171516E+0005

1000 5.63017732261513E+0005 6000 5.52669072420557E+0005

1500 6.48707174319903E+0005 6500 5.49492268154213E+0005

2000 7.31033980110939E+0005 7000 5.89081392532593E+0005

2500 9.85406499440432E+0005 7500 1.00681270240579E+0006

3000 1.13471345409721E+0006 8000 4.98314139143450E+0005

3500 4.95840415399220E+0005 8500 8.27839352551884E+0005

4000 4.22877980706847E+0005 9000 5.90302517716381E+0005

4500 5.32958603946419E+0005 9500 6.34971039929613E+0005

5000 6.12086751666785E+0005 10000 6.49020448431187E+0005

Значения нормы совершают ограниченные колебания, что свидетельствует об устойчивости.

) and )

) and ) and

; k:=0; end;

С помощью аттрактора Лоренца можно моделировать многие другие явления. Например, Хакен получил уравнения Лоренца, решая задачу нерегулярного распределения максимумов лазерного излучения, Йорке и Йорке обнаружили его, решая задачу о конвекции в тороидальной области. Кноблох нашел, что к системе Лоренца сводится задача о дисковом динамо. Педлоски и Френтен использовали уравнения Лоренца для описания динамики слабонеустойчивых барокли-нических волн конечной амплитуды. Существуют и другие задачи, которые можно моделировать этими уравнениями [3]. Систему Лоренца, по крайней мере на малых временах, можно использовать для описания динамики небольших городских систем, входящих в состав метрополии.

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

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

Предполагается, что локационные характеристики городского пространства описываются следующими тремя переменными: X - продукция, производимая городской системой, Y - численность коренного населения, Z - земельная рента.

Продукция городской промышленности может идти на потребление населения или экспортироваться вовне. В этих условиях возможна следующая модель динамика города [3]

= а (аУ - а-,Х), йг 1 2 3

ау

• -= с1(с2X - с3У) - с4XZ,

аг

— = аху - а ^.

„ аг 1 2

Последняя модель является системой Лоренца, что еще раз подчеркивает актуальность ее исследования.

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

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

2. Ромм Я.Е. , С.Г. Буланов Компьютерный анализ устойчивости по Ляпунову систем линейных дифференци-

альных уравнений. - Таганрог: ТГПИ имени А.П. Чехова, 2012. - 148 с.

3. Занг, В.Б.Синергетическая экономика. Время и перемены в нелинейной экономической теории. - М.: Мир

1999. - 335 с.