авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

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


Pages:     | 1 | 2 || 4 |

«Министерство науки и образования РФ Федеральное Государственное Бюджетное Образовательное Учреждение Высшего Профессионального Образования Уральский ...»

-- [ Страница 3 ] --

2.3.2. Интегральное эволюционное уравнение параметрической экономической системы кейнсианского типа Рассмотрим получение интегрального эволюционного уравнения парамет рической экономической системы в рамках модели Кейнса с непрерывным време нем, находящейся под воздействием внешних возмущающих факторов.

Динамика такой системы описывается возмущнной задачей Коши для уравнения (1.4.10) с начальным условием (1.4.2):

dyt 1 ct0 yt ct yt f t, (2.3.7) dt yt0 y0, (2.3.8) где в соответствие с (2.1.5) 1 d k ct ct y t t t0 k yt.

k k 1 k! dt В соответствии с параграфом 2.1 решение задачи Коши (2.3.7), (2.3.8) для имеющего экономический смысл начального условия y0 0 записывается в виде y t Gt s cs y s f s ds 1 ct0 y0 Z t, (2.3.9) где функция Грина определяется через решение специальной задачи Коши dZ t 1 ct0 Z t 0, dt Z 0 и для рассматриваемого случая имеет вид Gt H t e1ct0 t. (2.3.10) Подстановка (2.3.10) в (2.3.9) приводит к следующему интегральному эволюци онному уравнению t yt e 1c t0 t s f s cs ys ds 1 ct0 y0e 1c t0 t. (2.3.11) t Вводя обозначения Gt, s e 1c t0 t s и t y0 t Gt, s f s ds, t уравнение (2.3.11) перепишем в виде t yt y0 t 1 ct0 y0e e 1c t0 t s cs ys ds. (2.3.12) 1c t0 t t Теперь на основании следствия из утверждения 2.2.1 можно утверждать, что интегральное эволюционное уравнение для изучаемой системы (2.3.12) эквива лентно задаче Коши (2.3.7), (2.3.8) и, следовательно, описывает эволюцию эконо мической системы (однопродуктовой экономики или малого предприятия), нахо дящейся под воздействием внешних возмущающих факторов, в модели Кейнса.

Решение уравнения (2.3.12) численными методами находится значительно проще, чем решение задачи Коши (2.1.1), (2.1.2) стандартными методами.

Уравнение Вольтерра (2.3.12) является частным случаем интегрального уравнения Фредгольма. Действительно, рассматривая эволюцию системы на про межутке времени t0, T, видим, что уравнение (2.3.12) по своему смыслу опреде лено в нижней половине квадрата t0 t T, t0 s t декартовой плоскости. До определим ядро уравнения (2.3.12) тождественным нулм в верхней половине квадрата t0 t T, s t, положив Gt, s 0 для значений s t, то есть поло жим Gt, s t s e 1ct0 t s, где 1, s t;

t s 0;

s t.

Тогда верхний предел интегрирования можно положить равным T и записать уравнение (2.3.12) в квадрате t0 t T, t0 s T декартовой плоскости в виде уравнения Фредгольма второго рода:

T yt y0 t 1 ct0 y0e t s e 1c t0 t s cs ys ds.

1c t0 t t (2.3.13) Интегральное эволюционное уравнение (2.3.13), описывающее эволюцию одномерной экономической системы (экономики) во времени, является аналогом уравнения Липмана-Швингера (УЛШ) теории рассеяния, записанным, однако во временном представлении. По причине указанной аналогии функция ct (флук туация функции склонности к потреблению) в уравнении (2.3.13) может быть на звана “потенциалом склонности к потреблению”.

Отметим частный случай прямой динамической задачи, когда склонность к t t0, T pt в уравнении (2.1) постоянна, то есть накоплению pt p const. Легко видеть, что для этого случая решение задачи Коши dyt pyt f t, yt0 y0, dt датся формулой yt y0 t 1 ct 0 y0 e 1c t0 t, которая может быть получена, конечно, путм аналитического решения задачи Коши методом вариации произвольной постоянной [36, 39].

2.3.3. Интегральное эволюционное уравнение экономической системы в модели Самуэльсона-Хикса Рассмотрим получение интегрального эволюционного уравнения парамет рической экономической системы в рамках модели Самуэльсона-Хикса с непре рывным временем, находящейся под воздействием внешних возмущающих фак торов.

Динамика такой системы описывается возмущнной задачей Коши (1.3.17), (1.3.18):

d 2 y t dyt dyt 1 r t0 1 ct0 yt r t qt yt f t, dt dt dt (2.3.30) dyt yt0 y0, y1, (2.3.31) dt где флуктуации коэффициентов 1 d k pt0 1 d k r t pt pt pt 0 t t0 t t0 k, k k k k 1 k! k 1 k!

dt dt (2.3.32) 1 d k qt0 1 d k ct qt qt qt 0 t t0 t t0 k, k k k k 1 k! k 1 k!

dt dt (2.3.33) определяются через флуктуации параметров системы 1 d k r t r t t t0 k, (2.3.34) k k 1 k! dt 1 d kc ct t0 t t0 k. (2.3.35) k k 1 k! dt Применяя формулу (2.1.17), получаем dZ t t yt Z t s f s ds y0 p10 y0 y1 Z t. (2.3.36) dt t Подстановка всех входящих в формулу (2.3.36) величин приводит к интегрально му эволюционному уравнению следующего вида:

dZ t 1 r t0 y0 y1 Z t yt y0 t dt dys G t, s r s cs y s f s ds.

ds (2.3.37) Учитывая, что Gt, s H t s Z t, уравнение (2.3.37) перепишем в окончатель ном виде:

dZ t 1 r t0 y0 y1 Z t yt y0 t dt dys t Z t, s r s cs ys f s ds. (2.3.22) ds t Полученное интегральное уравнение (2.3.22) является уравнением Вольтер ра второго рода. Из приведнного вывода следует, что уравнение (2.3.22) эквива лентно задаче Коши (1.3.17), (1.3.18) или (2.3.30), (2.3.31). Поэтому динамику экономической системы можно изучать, решая интегральное эволюционное урав нение (2.3.22). Численное решение уравнения (2.3.22) достаточно просто можно реализовать методом последовательных приближений, сходимость которого га рантирована соответствующей теоремой для интегральных уравнений Вольтерра.

Так же, как уравнение (2.3.12), уравнение (2.3.22) является частным случаем интегрального уравнения Фредгольма. Действительно, рассматривая эволюцию системы на промежутке времени t0, T, видим, что уравнение (2.3.22) по своему смыслу определено в нижней половине квадрата t0 t T, t0 s t декартовой плоскости. Доопределим ядро уравнения (2.3.22) тождественным нулм в верхней половине квадрата t0 t T, s t, положив Gt, s 0 для значений s t, то есть положим Gt, s t s Z t, где 1, s t;

t s 0;

s t.

Тогда верхний предел интегрирования можно положить равным T и записать уравнение (2.3.22) в квадрате t0 t T, t0 s T декартовой плоскости в виде уравнения Фредгольма второго рода:

dGt 1 r t0 y0 y1 Gt yt y0 t dt dys T Gt, s r s cs ys f s ds, (2.3.23) ds t где использовано оговоренное выше соглашение об обозначении функции Грина.

Функция Грина линейного обыкновенного дифференциального оператора с постоянными коэффициентами 1 r t0 и 1 ct0 является решением задачи Кош d 2 Z t dZ t 1 r t0 1 ct0 Z t 0, (2.3.24) dt dt dZ Z 0 0, 1, (2.3.25) dt В зависимости от соотношения коэффициентов 1 r t0 и 1 ct0 дискриминант D 1 r t0 41 ct0 характеристического уравнения k 2 1 r t0 k 1 r t0 0 (2.3.26) будет строго положительным, строго отрицательным или равным нулю.

Если корни характеристического уравнения вещественные и различные, что имеет место при условии D 0, то общее решение уравнения (2.3.24) имеет вид 1r t0 1r t0 2 41c t0 1r t0 1r t0 2 41c t t t G(t ) C1e C2 e 2.

Подстановка в начальные условия (2.3.25) дат следующие значения произволь ных постоянных:

1 C1 ;

C2.

1 r t0 41 ct0 1 r t0 41 ct 2 Следовательно, функция Грина уравнения (2.13) в этом случае имеет вид:

1 r t0 1 r t0 2 41c t0 1 r t0 1 r t0 2 41c t t t e 2 e G (t ). (2.3.27) 1 r t0 41 ct Если корни характеристического уравнения вещественные и равные (харак теристическое уравнение имеет вещественный корень кратности 2), что имеет ме сто при условии D 0, то общее решение уравнения (2.3.24) имеет вид 1r t0 1r t 1 r t0 Gt C t t e C2 t e 2 2.

Подстановка в начальные условия (2.3.25) дат для произвольных постоянных следующие значения: C1 0 ;

C 2 1. Следовательно, функция Грина уравнения (2.3.24) в этом случае имеет вид:

1 r t t G (t ) te. (2.3.28) Если корни характеристического уравнения комплексно-сопряжнные, что имеет место при условии D 0, то общее решение уравнения (2.3.24) имеет вид 41 ct 1 r t 1 r t t C1 cos t G (t ) e 0 41 ct 1 r t C2 sin t.

0 Подстановка в начальные условия (2.3.25) дат следующие значения произволь ных постоянных: C1 0 ;

C2. Следовательно, функция 41 ct0 1 r t Грина уравнения (2.3.24) в этом случае имеет вид 41 ct 1 r t 1r t t sin t.

G (t ) e 0 41 ct0 1 r t 2 (2.3.29) Подставляя полученные выражения для функции Грина (2.3.27), (2.3.28), (2.3.29) последовательно в уравнение (2.3.23), получим интегральные эволюцион ные уравнения для трх случаев существования корней характеристического уравнения (2.3.26). В силу очевидного характера проводимой подстановки приво дить эти уравнения не будем.

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

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

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

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

2.3.4. Система интегральных эволюционных уравнений многомерной параметрической экономической системы в межотраслевой балансовой модели Леонтьева В случае, когда структурные параметры многомерной экономической сис темы можно считать гладкими функциями, вместо представления по формуле Тейлора, для них можно записать разложение в ряд Тейлора и дифференциальное уравнение (1.3.33) принимает вид:

d yt At yt Pt yt f t, I (2.3.30) dt где d k Pt0 t t k def Pt A Pt. (2.3.31) dt k k!

k Тогда эволюционное интегральное уравнение, эквивалентное возмущнной задаче Коши d yt At yt Pt yt f t, I (2.3.32) dt yt0 y0, (2.3.33) в соответствие с утверждением 2.2.1 для рассматриваемого случая примет вид t yt Gt, t0 y0 y0 t Gt, s Ps ys ds, (2.3.34) t t где вектор y0 t dsGt, s f s.

t Векторное интегральное эволюционное уравнение (2.3.34) можно перепи сать в виде системы уравнений Вольтерра y i t a k y 0 y 0 t Gk t, q k y j d, tn n n ik i i (2.3.35) a j 1 k 1 j k def t n y t Gk t, f k d, i i (2.3.36) a k где t a, b, t, а q ij t a ij p ij t.

Доопределяя матричные элементы Gk t, s в верхней половине t квад i рата a t, b условием i, k 1, n Gk t, 0, перепишем систему уравне i ний Вольтерра (2.3.35) для динамической межотраслевой балансовой модели в виде системы уравнений Фредгольма второго рода cn y t a y y t K ij t, y j d, n i i k i (2.3.37) k 0 k 1 a j где c a, b, а ядра K ij t, Gk t, q k.

n i j k ГЛАВА 3. ОБРАТНАЯ ЗАДАЧА ДЛЯ ПАРАМЕТРИЧЕСКОЙ СИСТЕМЫ 3.1. Постановка обратной задачи для параметрической системы 3.1.1. Задача определения коэффициентов системы дифференциальных уравнений Рассмотрим постановку и метод решения обратной задачи сразу для систе мы обыкновенных дифференциальных уравнений вида (1.3.3) 1 0 0 y1 t p1 t p1 t p1 t y1 t f 1 t 2 2 2 2 n 0 1 0 d y t p1 t p2 t pn t y t f t 2.......... dt....................,....................

............

0 0 1 y n t p n t p n t p n t y n t f n t 1 2 n (3.1.1) описывающей динамику параметрической системы с сосредоточенными парамет рами одного из предметных типов, рассмотренных в данной работе выше. Дальше систему уравнений (3.1.1) будем рассматривать как одно векторное обыкновенное дифференциальное уравнение, для которого и будет поставлена задача Коши.

Итак, пусть эволюция динамической системы описывается решением задачи Ко ши d y Pt y f, (3.1.2) dt yt0 y0. (3.1.3) Раскладывая матрицу Pt в ряд Тейлора в окрестности t 0, приведм систему уравнений (3.1.2) к виду (1.3.24) d yt A yt Pt yt f t, I (3.1.4) dt где согласно (1.3.25) флуктуационная матрица d k Pt0 t t k Pt A Pt. (3.1.5) dt k k!

k Векторное дифференциальное уравнение (3.1.4) описывает эволюцию во времени многомерной параметрической системы с сосредоточенными параметра ми, имеющей n входов и n выходов (рисунок 3.1.1), под воздействием внешних возмущений.

Входы системы Выходы системы f t yt S Рис. 3.1.1. Иллюстрация к постановке обратной динамической задачи для многомерной системы Во введении сформулирована обратная динамическая задача (ОДЗ) для ди намической системы, как задача нахождения оператора системы S при известном входном воздействии и известном отклике системы. В соответствие с обозначе ниями рисунка 3.1.1 обратная динамическая задача для многомерной параметри ческой системы может быть сформулирована так:

по известным на всм промежутке времени эволюции входным воздействи ям на систему f t, отклику системы yt и дифференциальному оператору фоновой системы найти дифференциальный оператор p1 t p1 t p1 t 1 0 0 2 2 n 0 1 0 d p1 t p2 t pn t 2 Lt.......... dt.................... (3.1.6)....................

............

0 0 1 p n t p n t p n t 1 2 n системы эволюционных уравнений (3.1.1).

Так сформулированная обратная динамическая задача является задачей оп ределения коэффициентов, точнее вариаций p ij t a ij p ij t коэффициентов системы дифференциальных уравнений (3.1.1). Напомним, что сформулированная постановка обратной задачи не исчерпывает всего множества возможных поста новок обратных задач. Второй, представляющей интерес для практики обратной задачей, является задача определения правой части уравнения (3.1.2) при извест ных элементах матрицы Pt и входном воздействии f t. Эта постановка об ратной задачи в работе не рассматривается.

3.1.2. Предметный смысл обратной динамической задачи Обсудим предметный (физический) смысл обратной динамической задачи определения коэффициентов дифференциального уравнения (3.1.1) или (3.1.2).

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

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

задавая входные воздействия и отклик системы на них, а определять отклонения параметров технической системы от фиксированных значений для фоновой сис темы.

Известно [31], что решение обратной задачи в случае, когда прямая задача формулируется как краевая (или граничная задача), может быть сведено к реше нию нелинейного интегрального уравнения первого рода или системы таких уравнений. Смысл такого “сведения” состоит в замене задачи определения диф ференциального оператора L, задачей определения интегрального оператора сис темы S, то есть обратной задачей в изначальной формулировке.

Как будет показано ниже, не смотря на то, что в нашем случае прямая зада ча является задачей Коши, обратная динамическая задача также может быть све дена к решению нелинейного интегрального уравнения первого рода, отличающе гося от уравнения для краевых задач, выписанного в работе [16]. Более того, на основе результатов работ [46, 47], посвящнных стационарной теории рассеяния волновых полей, ниже показано, что решение обратной динамической задачи для параметрической системы с сосредоточенными параметрами может быть получе но как решение некоторого линеаризованного интегрального уравнения, которое в упомянутых работах названо линеаризованным уравнением реконструкции.

3.2. Метод решения обратной динамической задачи для параметрической системы с сосредоточенными параметрами 3.2.1. Нелинейное интегральное уравнение обратной динамической задачи Рассмотрим теперь “интегральную” формулировку обратной динамической задачи для параметрической системы. В соответствие с обозначениями рисунка 3.1.1 эта формулировка принимает следующий вид:

по известным на всм промежутке времени эволюции входным воздействи ям на систему f t и отклику системы yt найти S – матрицу (оператор) сис темы S t.

Считая известной матричную функцию Грина для начальных условий (3.1.3), запишем интегральное эволюционное уравнение системы в виде (2.1.69) t yt Gt, t0 y0 y0 t Gt, s Ps ys ds, (3.2.1) t где вектор t y0 t dsGt, s f s. (3.2.2) t Дальше, не умаляя общности, рассмотрим задачу Коши с нулевыми началь ными условиями yt0 0. (3.2.3) Тогда интегральное эволюционное уравнение (3.2.1) примет несколько более про стой вид t yt y0 t Gt, s Ps ys ds. (3.2.4) t Интегральное уравнение (3.2.4) описывает эволюцию во времени парамет рической системы, имеющей n входов и n выходов (рисунок 3.1.1). Перепишем уравнение (3.2.4) в виде уравнения Фредгольма, вводя обозначения параграфа 1. b y t y0 t G t, s Ps y s ds, (3.2.5) a 0, s t;

где наложено условие Gk t, s i Внеинтегральный член (3.2.2) в соот 0, s t.

ветствии с этим условием имеет вид b y0 t A y0 dsGt, s f s.

a Основываясь на результатах параграфа 1.4, запишем решение интегрально го эволюционного уравнения (3.2.4) (решение задачи Коши (3.1.2), (3.1.3)) в виде борновского ряда b y t y0 t G t, t1 Pt1 y0 t1 dt a bb G t, t1 Qt1 G t1, t 2 Pt 2 y0 t 2 dt2 dt aa bbb G t, t1 Pt1 G t1, t 2 Pt 2 G t 2, t3 Pt3 y0 t3 dt3 dt2 dt1.

aaa (3.2.6) Как было показано в главе 1, эволюция системы может быть выражена явно через е S – матрицу (оператор системы) по формуле (1.4.4) yt S y0 t, (3.2.7) где матрица S I T, а матрица (оператор) взаимодействия b bb T dt1G t, t1 Pt1 dt1dt2G t, t 2 Pt 2 G t 2, t1 Pt a aa bbb dt1dt2 dt3G t, t3 Pt3 G t3, t 2 Pt 2 G t 2, t1 Pt1.

aaa Уравнение (3.2.7) запишем в виде y t, p ij S t, p ij y0 t, (3.2.8) где указана явная зависимость S – матрицы и отклика системы от элементов мат рицы вариаций параметров системы.

Для решения обратной задачи уравнение (3.2.8) нужно представить в дуаль ной форме y t, p ij S t, y0 t p ij, (3.2.9) где p ij – управляющие параметры системы. Тогда решение обратной задачи сво дится к обращению оператора S и записывается в виде p ij S 1 y. (3.2.10) 3.2.2. Интегральные уравнения реконструкции Для обращения оператора системы применим метод Ньютона. В результате обращения оператора системы будет получена система интегральных уравнений реконструкции, решение которой при известных входному воздействию на систе му и отклику на него позволит найти значения управляющих параметров, обеспе чивающие достижение требуемого отклика.

Вводя обозначение def y t, p ij y0 t, S i y t, p (3.2.10) j запишем борновский ряд (3.2.6) для рассматриваемой системы в виде y t, p b dt1G t, t1 Pt S i j a bb dt1dt2G t, t 2 Pt 2 G t 2, t1 Pt aa bbb dt1dt2 dt3G t, t3 Pt3 G t3, t 2 Pt 2 G t 2, t1 Pt1 y0 t1.

aaa (3.2.11) Здесь скобкой обозначено место, на которое в интегральном операторе по мещается вектор-функция y0 t1. Если считать входное воздействие y0 t1 на систему известным, а отклик системы t S t, p ij задать (или измерить экспери ментально), то уравнение (3.2.11) – это нелинейное интегральное уравнение пер вого рода типа Ляпунова-Шмидта, записанное в дуальной форме [5], то есть, от носительно неизвестных управляющих параметров p ij t i, j 1, 2, 3,, n.

Очевидно, что решение этого уравнения можно найти только численными методами. Имеющееся аналитическое представление решения подобных уравне ний в виде ряда последовательных приближений [41] не представляет интереса в виду его чрезмерной громоздкости. В силу специальной структуры уравнения (3.2.11) весьма привлекательно, на первый взгляд, выглядит попытка его решения методом дискретизации, в результате которого векторное интегральное уравнение (3.2.11) сведтся к системе нелинейных алгебраических уравнений. Решение сис темы нелинейных алгебраических уравнений, однако, является весьма сложной и громоздкой задачей и по этой причине этот путь решения уравнения (3.2.11) явля ется, по-видимому, тупиковым.

В работах [46, 47] показано, что путм линеаризации векторное нелинейное интегральное уравнение (3.2.11) можно свести к системе линейных интегральных уравнений первого рода стандартного вида, которые названы линеаризованными уравнениями реконструкции (ЛУР), так как посредством их решения можно осу ществить восстановление управляющих параметров (реконструкцию) исследуе мого объекта. Приведм упрощнный вывод ЛУР для исследуемого случая ли нейной многомерной системы с сосредоточенными параметрами, имея в виду применение теории к соответствующей модели экономической или радиотехни ческой системы.

Перепишем уравнение (3.2.11) в усечнном виде, ограничившись на время процессами взаимодействия третьей кратности:

y S t, p ij dt1Gt, t1 Pt1 y0 t b a bb dt1dt2G t, t 2 Pt 2 Gt 2, t1 Pt1 y0 t aa bbb dt1dt2 dt3G t, t3 Pt3 G t3, t 2 Pt 2 G t 2, t1 Pt1 0 t1.

aaa (3.2.12) y 0 t, где t a, b, системы на некоторое S Допустим, что известен отклик управляющее воздействие P0 t нулевого приближения. Тогда уравнение (3.2.12) для известного отклика и управляющих параметров запишется в виде b y t, P0 t dt1G t, t1 P0 t1 y0 t S a bb dt1dt2Gt, t 2 P0 t 2 G t 2, t1 P0 t1 0 t aa bbb dt1dt2 dt3G t, t3 P0 t3 G t3, t 2 P0 t 2 G t 2, t1 P0 t1 0 t1.

aaa (3.2.13) Уравнение (3.2.12), вводя обозначение Ptk P tk Ptk, где Ptk – ва риации управляющих параметров Ptk, перепишем в виде b y t, P dt1Gt, t1 P0 t1 Pt1 y0 t S a bb dt1dt2G t, t 2 P0 t 2 Pt 2 G t 2, t1 P0 t1 Pt1 y0 t aa bbb dt1dt2 dt3G t, t3 P0 t3 Pt3 G t3, t 2 P0 t 2 Pt aaa Gt2, t1 P0 t1 Pt1 y0 t1.

(3.2.14) Вычтем уравнение (3.2.13) из уравнения (3.2.14). Обозначая y S t, P y S t, P y0 t, P0, S запишем результат вычитания:

y S t, P b dt1G t, t1 P t1 0 t a dt1 dt2 G t, t 2 P0 t 2 G t 2, t1 Pt1 G t, t 2 Pt 2 G t 2, t1 P0 t bb aa Gt, t 2 Pt 2 Gt 2, t1 Pt1 y0 t dt1 dt2 dt3 G t, t 3 P0 t 3 GGt 3, t 2 P0 t 2 G t 2, t1 P t bbb aaa G t, t 3 P0 t 3 G t 3, t 2 Pt 2 G t 2, t1 P0 t1 G t, t 3 P0 t 3 G t 3, t 2 Pt 2 G t 2, t1 Pt G t, t 3 Pt 3 G t 3, t 2 P0 t 2 G t 2, t1 P0 t1 G t, t 3 Pt 3 G t 3, t 2 P0 t 2 G t 2, t1 Pt G t, t 3 Pt 3 G t 3, t 2 Pt 2 G t 2, t1 P0 t Gt, t3 Pt3 Gt3, t 2 Pt 2 Gt 2, t1 Pt1 y0 t1. (3.2.15) В правой части уравнения (3.2.15) первое слагаемое соответствует однократному, второе – двукратному, а третье – трхкратному взаимодействию. Уравнение (3.2.15) – нелинейное интегральное уравнение первого рода. Проведм его линеа ризацию. Для этого выделим в уравнении (3.2.15) члены, начиная со второго по рядка малости по вариациям P t :

y S t, P b dt1G t, t1 P t1 y 0 t a dt1 dt2 G t, t 2 P0 t 2 G t 2, t1 Pt1 G t, t 2 Pt 2 G t 2, t1 P0 t1 y 0 t bb aa dt1 dt2 dt3 G t, t 3 P0 t 3 G t 3, t 2 P0 t 2 G t 2, t1 Pt bbb aaa G t, t 3 P0 t 3 G t 3, t 2 Pt 2 G t 2, t1 P0 t Gt, t3 Pt3 Gt3, t 2 P0 t 2 Gt 2, t1 P0 t1 y0 t члены, начиная со второго порядка по P. (3.2.16) В соответствии с методом Ньютона пренебрежм в уравнении (3.2.16) чле нами, начиная со второго порядка по вариациям P t. В результате такого пре небрежения получим искомое линеаризованное уравнение реконструкции:

y S t, P b dt1G t, t1 P t1 y 0 t a dt1 dt2 G t, t 2 P0 t 2 G t 2, t1 Pt1 G t, t 2 Pt 2 G t 2, t1 P0 t1 y 0 t bb aa dt1 dt2 dt3 G t, t 3 P0 t 3 G t 3, t 2 P0 t 2 G t 2, t1 P t bbb aaa G t, t 3 P0 t 3 G t 3, t 2 Pt 2 G t 2, t1 P0 t Gt, t3 Pt3 Gt3, t 2 P0 t 2 Gt 2, t1 P0 t1 y0 t1. (3.2.17) Полученное уравнение является приближнным. Приближнный характер уравнения (3.2.17) обусловлен, с одной стороны, кратностью учитываемых про цессов взаимодействия в системе, а с другой стороны – величиной отброшенных членов высших порядков по вариациям P t. Для того, чтобы ЛУР (3.2.17) можно было эффективно решать, его следует преобразовать к стандартному для линей ных интегральных уравнений виду. Обозначим входящие в ЛУР (3.2.17) интегра лы так:

b I11 dt1G t, t1 Pt1 y0 t1 ;

(3.2.18) a bb I 21 dt1dt2G t, t 2 P0 t 2 G t 2, t1 Pt1 y0 t1 ;

(3.2.19) aa bb I 22 dt1dt2G t, t 2 Pt 2 G t 2, t1 P0 t1 y0 t1 ;

(3.2.20) aa bbb I 31 dt1dt2 dt3Gt, t3 P0 t3 Gt3, t 2 P0 t 2 Gt 2, t1 Pt1 y0 t1 ;

aaa (3.2.21) bbb I 32 dt1dt2 dt3Gt, t3 P0 t3 Gt3, t 2 Pt 2 Gt 2, t1 P0 t1 y0 t1 ;

aaa (3.2.22) bbb I 33 dt1dt2 dt3Gt, t3 Pt3 Gt3, t 2 P0 t 2 Gt 2, t1 P0 t1 y0 t1.

aaa (3.2.23) Теперь ЛУР (3.2.17) переписывается в символическом виде S t I 11 I 21 I 22 I 31 I 32 I 33. (3.2.24) Интегралы в уравнении (3.2.24) преобразуем последовательно.

Интеграл (3.2.18), осуществляющий учт процессов однократного взаимо действия, не преобразуется.

Интегралы (3.2.19) и (3.2.20), осуществляющие учт процессов двукратного взаимодействия преобразуются, соответственно, к виду bb I 21 dt1dt2Gt, t 2 P0 t 2 G t 2, t1 Pt1 y0 t aa b dt1G 1 t, t1 Pt1 y0 t1, (3.2.25) a bb I 22 dt1dt2G t, t 2 Pt 2 G t 2, t1 P0 t1 y0 t aa b dt1G t, t1 Pt1 y1 t1. (3.2.26) a Здесь введены следующие обозначения:

def b t, t1 dt2Gt, t2 P0 t2 Gt2, t G (3.2.27) a – оператор распространения (пропогатор, или передаточная функция) первого по рядка;

def b y1 t1 dt2G t1, t 2 P0 t 2 y0 t 2 (3.2.28) a – входное воздействие, подвергнутое процессу однократного взаимодействия. Для краткости будем называть член вида (3.2.28) однократным взаимодействием.

Интегралы (3.2.21), (3.2.22) и (3.2.23), осуществляющие учт процессов трхкратного взаимодействия, преобразуются, соответственно, к виду bbb I 31 dt1dt2 dt3G t, t3 P0 t3 G t3, t 2 P0 t 2 G t 2, t1 Pt1 y0 t aaa b dt1G 2 t, t1 Pt1 y0 t1, (3.2.29) a bbb I 32 dt1dt2 dt3G t, t3 P0 t3 G t3, t 2 Pt 2 G t 2, t1 P0 t1 y0 t aaa b dt1G 1 t, t1 Pt1 y1 t1, (3.2.30) a bbb I 33 dt1dt2 dt3G t, t3 Pt3 G t3, t 2 P0 t 2 G t 2, t1 P0 t1 y0 t aaa b dt1G t, t1 Pt1 y2 t1. (3.2.31) a Здесь введены обозначения def b t, t1 dt2G 1 t, t2 P0 t2 Gt2, t1, G (3.2.32) a def b y2 t1 dt2G t1, t 2 P0 t 2 y1 t 2, (3.2.33) a для пропогаторов второго порядка и двукратных взаимодействий, соответственно.

Подставляя (3.2.25) – (3.2.33) в уравнение (3.2.24), после группировки инте гралов получаем:

y t dt1 0 t, t1 Pt1 y0 t1 y1 t1 y2 t1 b S G a G 1 t, t1 Pt1 y0 t1 y1 t1 G 2 t, t1 Pt1 y0 t1.

(3.2.34) Здесь для единообразия введено обозначение G 0 t, t1 Gt, t1. Легко видеть, что уравнение (3.2.34) с учтом известного свойства сопряжнных матриц At yt yt A t, можно переписать в следующем виде:

b y t dt1 G 0 t, t1 y 0 t1 y 1 t1 y 2 t1 P t S a G 1 t, t1 y 0 t1 y 1 t1 P t1 G 2 t, t1 y 0 t1 P t1.

(3.2.35) В уравнении (3.2.35) вектор yt y1 t y 2 t y n t – транспонирован ный по отношению к вектору-столбцу вектор-строка.

Уравнение (3.2.35) является искомым ЛУР, учитывающим процессы взаи модействия в системе до третьей кратности включительно. Все выполненные при выводе уравнения (3.2.35) преобразования справедливы для любой конечной кратности взаимодействия. Поэтому ЛУР (3.2.35) легко обобщается и для лю бой конечной кратности взаимодействия принимает вид:

1 b y t, P dt1 G t, t1 y t1 P t1, S (3.2.36) 0 a где пропогаторы и кратные взаимодействия различных порядков рассчитываются по формулам b t, t1 dt' G 1 t, t 'P0 t 'G 0 t ', t1, 1 1, G a b y t1 dt ' G 0 t1, t 'P0 t ' y 1 t '.

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

i j Итерации продолжаются, пока невязка левой части не достигнет заданного значе ния.

Отметим, что ЛУР (3.2.36) записано в операторной форме. В матричной форме записи ЛУР (3.2.36) эквивалентно системе линейных скалярных инте гральных уравнений первого рода. Задача решения таких уравнений является примером некорректно поставленной задачи. Для построения устойчивого алго ритма решения ЛУР (3.2.36) следует применять один из классических методов ре гуляризации, например, метод М. М. Лаврентьева [31].

3.2.3. Замечания о природе приближнного характера ЛУР При выводе ЛУР в предыдущем пункте использована процедура линеариза ции нелинейного интегрального уравнения (3.2.11). Поясним смысл этой проце дуры. Уравнение (3.2.11) можно переписать в дуальной форме y S t, P S t, P Pt, (3.2.37) где оператор системы S t, P имеет вид b S t, P I dt1G t, t1 Pt a bb dt1dt2G t, t 2 Pt 2 G t 2, t1 Pt aa bbb dt1dt2 dt3G t, t3 Pt3 G t3, t 2 Pt 2 G t 2, t1 Pt1.

aaa По определению оператор S t, P дифференцируем в некоторой окрестности точки P0, если в этой окрестности справедливо представление y S t, P y S t, P0 Dt, P0 Pt P, (3.2.38) где Dt, P – ограниченный линейный оператор, а P – бесконечно ма лые члены, начиная со второго порядка по совокупности управляющих парамет ров Pt. Линейный оператор Dt, P называется дифференциалом Фреше [66] оператора S t, P в дуальной записи (3.2.37). Линеаризация уравнения (3.2.11) означает пренебрежение в уравнении членами, начиная со второго поряд ка малости по вариациям P, то есть членами P. Таким образом, урав нение (3.2.36) является существенно приближнным уравнением, учитывающим, однако, процессы многократных взаимодействий в системе.

ГЛАВА 4. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ЛИНЕЙНОЙ ПАРАМЕТРИЧЕСКОЙ СИСТЕМЫ 4.1. Программная реализация алгоритма численного моделирования динамики линейной параметрической системы 4.

1.1. Структура программного комплекса для численного моделирования эволюции параметрической системы Программный комплекс для численного моделирования динамики линейной параметрической системы с сосредоточенными параметрами (“PARSIS”) пред ставляет собой унифицированную вычислительную программу для численного решения первой основной задачи, основанную на векторном представлении мате матической модели многомерной параметрической системы (2.1.68), (2.1.69), по лученном с помощью матричной функции Грина (2.1.55), доказанных в главе утверждениях 2.1.1, 2.1.2, 2.2.1 и следствия из утверждения 2.2.1 и конструкции S -матрицы системы (2.2.2) – (2.2.5). Программа является унифицированной, так как предназначена для математического моделирования как одномерных, так и многомерных параметрических систем с сосредоточенными параметрами. Про грамма PARSIS написана на языке высокого уровня FORTRAN PS с использова нием лицензионного экземпляра программного продукта “FORTRAN 10.1”. С точки зрения структурного программирования программный комплекс представ ляет собой совокупность подпрограмм-процедур, объединнных в одно целое ба зовой (управляющей) программой, реализующей управление процессом ввода данных, численного решения системы интегральных уравнений Вольтерра второ го рода, вывода результатов моделирования “на печать”. Численное решение про изводится методом последовательных приближений, полностью эквивалентным методу подстановок, сходимость которого для системы интегральных уравнений Вольтерра гарантирована теоремой 2.2.1. Структурная схема программного ком плекса показана на рисунке 4.1.1.

Приведм краткое описание структуры программного комплекса и проведм краткий анализ его работы с позиций системного анализа.

Для проведения моделирования динамики конкретной параметрической системы в управляющую программу из файла «start.dat» заносятся следующие па раметры:

– размерность системы (идентификатор “nsis”);

– размерность разбиения промежутка интегрирования (идентификатор “nksi”);

– размерность разбиения промежутка времени эволюции системы (иденти фикатор “ntau”);

– значение безразмерного параметра точности вычислений (идентификатор EPS);

– начальный (идентификатор LEFT) и конечный (идентификатор RFGHT) моменты времени;

– массив значений элементов фоновой матрицы коэффициентов (идентифи катор FOM).

На вход основной программы податся поток входных данных – массив компонентов вектора входного сигнала (идентификатор XKSI). Компоненты этого вектора в зависимости от решаемой задачи могут иметь различный физический, экономический и т. д. смысл.

Размерности всех массивов данных задаются декларативным оператором REAL.

Управляющая программа заносит полученные данные в память компьютера и использует их в качестве параметров при вызове подпрограмм-процедур типа SUBROUTINE, посредством которых осуществляются все необходимые вычисле ния.

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

Далее управляющая программа организует итерационный процесс для ре шения векторного эволюционного интегрального уравнения (системы интеграль ных эволюционных) уравнений Вольтерра (2.1.68) t yt Gt, t0 y0 y0 t Gt, s Ps ys ds, (4.1.1) t эквивалентного задаче Коши (2.1.61) d yt Pt yt f t, yt0 y0.

I (4.1.2) dt Итерационный процесс реализуется посредством блока обратной связи и произ водится до тех пор, пока безразмерный параметр точности вычислений EPS не примет значение, меньшее заданного.

После окончания итерационного процесса программа записывает результа ты вычислений в файл «vent.dat».

SIGNAL FORCE KERNEL Вектор входного Вектор сигнала выходного сигнала Итерированный вектор входного Блок Управляющая сигнала обратной программа связи Декларации GRIN QVAR PSIS Рис. 4.1.1. Структурная схема программного комплекса Связь между управляющей программой PARSIS и подпрограммами SUB ROUTINE осуществляется посредством аппарата формальных параметров, в каче стве которых используются простые переменные и массивы. Для удобства ис пользуются массивы с переменными границами.

4.1.2. Подпрограммы-процедуры Приведм краткое описание назначения подпрограмм-процедур типа SUB ROUTINE.

1. Подпрограмма-процедура «CRASH» выполняет расчт узлов равномер ной сетки интегрирования и точек, в которых рассчитываются компоненты векто ра выходного сигнала.

При решении методом последовательных приближений векторного инте грального уравнения Вольтерра второго рода (4.1.1) значения верхнего предела фиксируются, то есть, производится разбиение промежутка интегрирования в двух уровнях. На рисунке 4.1.2 показано двухуровневое разбиение оси изменения переменной t в частном случае равномерного разбиения.

s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s A 0 1 3 Ось изменения переменной t..........

.............

Рис. 4.1.2. Разбиение оси изменения незави L симой переменной В каждом частичном промежутке разбиения k 1, k k 0, 1, 2, (рисунок 4.1.2) компоненты вектор-функции y s вычисляются путм решения векторно го интегрального уравнения Фредгольма второго рода, к которому сводится урав нение Вольтерра (4.1.1) при наложении на ядро уравнения условий вида [30] G s, если s t;

Gt s 0, если s t.

2. Подпрограмма-процедура «SUBROUTINE GRIN» выполняет расчт эле ментов матричной функции Грина (для начальных условий) для фоновой системы по формуле (2.1.55) Gt H t Z t H t Y t Y 1 t0. (4.1.3) 3. Подпрограмма-процедура «SUBROUTINE PSIS» выполняет расчт мат рицы коэффициентов Pt системы дифференциальных уравнений (4.1.2) для ка ждой точки разбиения промежутка времени, на котором исследуется эволюция системы.

4. Подпрограмма-процедура «SUBROUTINE MATVAR» выполняет расчт элементов матрицы вариаций параметров по формуле (2.1.67) P t A P t. (4.1.4) 4. Подпрограмма-процедура «SUBROUTINE FORCE» выполняет расчт компонентов вектора силового воздействия на систему.

5. Подпрограмма-процедура «SUBROUTINE SIGNAL» выполняет расчт компонентов вектора входного сигнала по формуле b y0 t G t, t0 y0 dsGt, s f s. (4.1.5) a 6. Подпрограмма-процедура «SUBROUTINE KERNEL» выполняет расчт матричного ядра векторного интегрального уравнения по формулам K t, s Gt, s Ps. (4.1.6) 7. Подпрограмма-процедура «SUBROUTINE SUMMA» выполняет числен ное интегрирование для каждого номера итерации.

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

Для математического моделирования одномерных параметрических систем с сосредоточенными параметрами разработаны также более простые вычисли тельные программы «ECHO» и программа «Direct», осуществляющие численное решение интегрального уравнения Вольтерра второго рода соответственно разно стным методом и методом последовательных приближений.

4.2. Примеры численного решения основной задачи для одномерной модели параметрической системы В этом параграфе приведены некоторые примеры численного моделирова ния (решения основной задачи) параметрических моделей одномерных экономи ческих систем. Подчеркнм, что по объективным причинам расчты носят чисто модельный характер. Результаты численного моделирования были опубликованы в работах [50 – 52].

4.2.1. Предварительные замечания об алгоритме численного решения основного интегрального уравнения Во второй главе работы показано, что решение возмущнной задачи Коши, описывающей эволюцию одномерной параметрической системы, находящейся под воздействием возмущающих факторов, может быть сведено к эволюционному интегральному уравнению, которое для системы, описываемой уравнением 1-го порядка или уравнением второго порядка, имеет соответственно вид (1.2.16) или (1.2.23). Во второй главе эти уравнения преобразованы в интегральные уравнения, описывающие динамику экономической системы в модели Кейнса (уравнение (2.3.12)) t yt y0 t 1 ct0 y0e e 1c t0 t s cs ys ds, 1c t0 t (4.2.1) t и в модели Самуэльсона-Хикса (уравнение (2.3.23)) dZ t 1 r t0 y0 y1 Z t yt y0 t dt dys t Z t, s r s cs ys f s ds. (4.2.2) ds t Толкование всех величин, входящих в уравнения (4.2.1) и (4.2.2), приведено в па раграфах 1.3 и 2.3. В параграфе 2.2 (теорема 2.2.1) доказана абсолютная и равно мерная сходимость ряда последовательных приближений (подстановок) для ре шения векторного интегрального эволюционного уравнения (2.2.1), частным слу чаем которого является как уравнение (4.2.1), так и уравнение (4.2.2).

Для численного решения одномерных интегральных уравнений, однако, можно применить и разностный метод. В работе [25] отмечено следующее: при менение разностного метода для решения одномерных интегральных уравнений приводит к хорошим результатам, а применение метода к двумерным (система из двух уравнений) интегральным уравнениям приводит к удовлетворительным ре зультатам. Так как тестирование программ численного решения интегральных эволюционных уравнений для моделей параметрических систем по объективным причинам затруднительно, то для такого тестирования применяется сравнение ре зультатов, полученных при использовании программ, реализующих различные алгоритмы.

В этом параграфе изложение разобьм на следующие два этапы.

1. На первом этапе будет проведено сравнение результатов численного ре шения интегрального уравнения (4.2.1) разностным методом и методом последо вательных приближений с точным аналитическим решением для модельной зада чи.

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

4.2.2. Разностный алгоритм решения интегрального эволюционного уравнения одномерной параметрической модели Рассмотрим техническую сторону простейшего алгоритма численного ре шения интегрального эволюционного уравнения [25].

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

Итак, рассматриваем интегральное эволюционное уравнение экономической сис темы в модели Кейнса следующего вида:

t yt y0 t Gt, p y d, (4.2.3) t где t y0 t Gt, f d. (4.2.4) t Здесь входной сигнал f t считается известной непрерывной функцией времени t. Подставляя (3.2) в (3.1), приведм интегральное уравнение к виду b y t G t, f p y d. (4.2.5) a В простейшем варианте конечно-разностного метода в уравнении (4.2.5) можно воспользоваться, например, непосредственно интегральной суммой Рима на (метод прямоугольников). Для этого, полагая t0 a, t b, рассмотрим разбие ние промежутка интегрирования a, b точками t j :

a, b a t0 t1 t2 tk 1 tk tn1 tn b.

В каждом частичном промежутке разбиения выберем среднюю точку, задав е формулой t j t j j 1, 2,, n.

x j t j 1 (4.2.6) Искомую функцию заменим ступенчатой функцией [15], то есть функцией, при нимающей постоянное значение на каждом частичном промежутке разбиения t, t j, положив j y j yx j j 1, 2,, n. (4.2.7) В уравнении (4.2.5) интеграл заменим интегральной суммой Римана y xi G xi, x j f x j px j y x j x j, n (4.2.8) j где x j t j t j 1 i, j 1, 2,, n. Получили систему линейных алгебраиче искомой функции в ских уравнений (СЛАУ) относительно значений y j y x j j 1, 2,, n разбиения промежутка времени эволюции системы.

точках x j Для эффективного решения СЛАУ (4.2.8) е нужно привести к стандартно му виду. Для упрощения записи введм дополнительные обозначения:

Gij Gxi, x j ;

f j f x j ;

j x j ;

p j px j.

Запишем для СЛАУ (4.2.8) следующую цепочку преобразований:

n n n n yi Gij f j j Gij p j j y j yi Gij p j j y j Gij f j j j 1 j 1 j 1 j n n n ij y j Gij p j j y j Gij f j j, j 1 j 1 j 1, i j, где ij Теперь СЛАУ (4.2.8) переписывается в виде:

0, i j.

y n n Gij f j j.

Gij p j j (4.2.9) ij j j 1 j Обозначим n ij ij Gij p j j, i Gij f j j. (4.2.10) j Подстановка (4.2.10) в (4.2.9) дат окончательный вид СЛАУ n ij j i i, j 1, 2,, n, (4.2.11) j или в разврнутом виде 111 122 1nn 1,, 21 1 22 2 2n n (4.2.12)................................................., n11 n 22 nnn n.

Система (4.2.12) имеет стандартный вид, и е решение может быть найдено при помощи одного из существующих численных методов, например, метода Га усса. Координаты вектора решений СЛАУ (4.2.12) y1, y2,, yn представляют T собой значения ВВП в различные моменты времени эволюции системы, принад лежащие множеству точек рассматриваемого разбиения промежутка времени эво люции.

4.2.3. Программная реализация разностного и итерационного алгоритмов решения интегрального эволюционного уравнения для одномерной параметрической модели Для численной реализации рассмотренного в предыдущем пункте алгорит ма решения интегрального эволюционного уравнения разработана вычислитель ная программа на языке Фортран, имеющая в соответствии с основными принци пами структурного программирования ярко выраженную структуру: “основная программа + подпрограммы”, реализующая разностный алгоритм решения инте грального уравнения. Структурная схема программы приведена на рисунке 4.2.1.

На вход основной программы подаются следующие величины:

1) “входной сигнал” – значения функции f t C I t, которую условно назо вм функцией инвестиций;

2) значение постоянного коэффициента p0 – склонность к накоплению для моде ли нулевого приближения;

3) размерность n сетки разбиения промежутка времени эволюции системы;

4) начальное и конечное значения времени эволюции системы.

Координаты точек сетки разбиения промежутка времени эволюции рассчи тываются в подпрограмме-процедуре “сетка” по формуле (4.2.6) и возвращаются в основную программу.

Значения функции Грина в точках сетки a, b a, b рассчитываются в подпрограмме-процедуре “Функция Грина” по формуле Gt, t1 H t aexp p0 t t1, приведнной в первой главе, и возвращаются в основную программу.

Значения коэффициентов СЛАУ рассчитываются в подпрограмме процедуре “Коэффициенты СЛАУ” по первой из формул (4.2.10) и возвращаются в основную программу.

Значения координат вектора правых частей СЛАУ (3.10) рассчитываются в подпрограмме-процедуре “Правые части СЛАУ” по второй из формул (4.2.10) и возвращаются в основную программу.

Метод Коэффици решения Сетка енты СЛАУ СЛАУ Входной Отклик ОСНОВНАЯ ПРОГРАММА сигнал системы Правые Функция части Грина СЛАУ Рис. 4.2.1. Структурная схема программы, реализующей разностный алгоритм Сформированная СЛАУ решается, например, методом Гаусса в подпро грамме “Метод решения СЛАУ”, после чего результаты расчтов возвращаются в основную программу и выводятся на печать.

Для чтения массивов входных данных использованы операторы форматного ввода. Связь между основной программой и подпрограммами-процедурами осу ществляется посредством передачи массивов формальных параметров, для кото рых использовано представление с переменными границами, что позволило сде лать программу более гибкой и удобной для моделирования. Текст основной про граммы приведн ниже.

Program Echo REAL del(5,5),cp(5), eou(5), ein(5), xp(5), alf(5,5),bet(5), gf(5,5) open(1, file="incom.dat") open(2, file="echo.dat") read(1,40) nd read(1,41) bl, br, p read(1,41) (cp(i), i=1,nd) read(1,41) (ein(i), i=1,nd) 40 format(i3) 41 format(5f10.3) do 1 i=1,nd do 1 j=1,nd del(i,j)=0.

if(i.eq.j)del(i,j)=1.

1 continue call break(nd, bl, br, xp,dx) call grin(nd, p0, xp, gf) call matrix(nd, del, gf, cp, dx, alf) call right(nd, gf, ein, dx, bet) call migc(nd, alf, bet) do 4 i=1,nd eou(i)=bet(i) 4 continue write(2,42) 42 format(5x"---ОТКЛИК СИСТЕМЫ---") write(2,41) (eou(i), i=1,nd) stop end Алгоритм последовательных приближений для решения интегрального эво люционного уравнения параметрической модели одномерной системы является частным случаем соответствующего алгоритма для многомерной модели, кото рый был подробно рассмотрен в первой главе работы. По этой причине мы его рассматривать не будем. Для численной реализации алгоритма последовательных приближений разработана вычислительная программа на языке Фортран, имею щая в соответствии с основными принципами структурного программирования ярко выраженную структуру: “основная программа + подпрограммы-процедуры”.

Структура программы повторяет основные черты программы численного модели рования многомерной системы, описанной в параграфе 6.1. Структурная схема программы приведена на рисунке 4.2.2.

Функция Сумма Сетка Грина Входной Отклик ОСНОВНАЯ ПРОГРАММА системы сигнал Итерационный цикл Рис. 4.2.2. Структурная схема программы, реализующей алгоритм последовательных приближений На вход основной программы податся следующие параметры задачи:

1. “входной сигнал” – значения функции f t C I t, которая условно на зывается функцией инвестиций;

2. значение постоянных коэффициентов p0 и q0 – коэффициента замедле ния и склонности к накоплению для фоновой модели;

3. размерность n сетки разбиения промежутка времени эволюции системы;

4. начальное и конечное значения времени эволюции системы.

Координаты точек сетки разбиения промежутка времени эволюции рассчи тываются в подпрограмме-процедуре “Сетка” по формуле (4.2.6) и возвращают ся в основную программу.

Значения функции Грина в точках сетки разбиения области a, b a, b в 0 зависимости от соотношения постоянных коэффициентов p1 и p 2 рассчитыва ются в подпрограмме-процедуре “Функция Грина” по формулам (2.2.16) – (2.2.18), полученным во второй главе работы, и возвращаются в основную про грамму.

В основной программе насчитываются значения ядра интегрального урав нения K ti, j G ti, j q j в точках сетки разбиения области a, b a, b.

На каждом шаге последовательных приближений интеграл рассчитывается по формуле I ti K ti, j y j x j в подпрограмме-процедуре “Сумма”.

n j Процесс последовательных приближений в программе реализован при помощи неявного цикла. В качестве опорного выбран входной сигнал y0 t, вычисляемый по формуле (4.2.4).

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

В силу доказанной в теореме 1.4.1 сходимости метода последовательных приближений, для окончания работы программы можно использовать простей ший критерий останова по величине евклидовой нормы разности значений функ ции отклика системы n-го и n 1 -го приближений.

Текст основной программы приведн ниже.

Program Direct dq(n) – управляющий параметр c eou(n) – отклик системы c ein(n) – входной сигнал c xp(n) – точки сетки разбиения c gf(n,n) – функция Грина c bl – левый конец промежутка времени эволюции c br - левый конец промежутка времени эволюции c с p1, p2 – коэффициенты уравнения для системы нулевого приближения fpr(n) – предыдущее приближение c fsu(n) – последующее приближение c REAL dq(10), f0(10), fou(10), fin(10), xp(10), gf(10,10) REAL ker(10,10), fpr(10), fsu(10) open(1, file="echo-in.dat") open(2, file="echo-out.dat") read(1,40) nd read(1,50) pogr read(1,41) bl, br, p1, p c read(1,41) (dq(i), i=1,nd) c read(1,41) (fin(i), i=1,nd) 40 format(i3) 41 format(5f10.3) 50 format(f10.5) call break(nd, bl, br, xp, dx) call signal(nd, xp, fin) call manage(nd, xp, dq) call grin(nd, p1, p2, xp, gf) c вычисление ядра интегрального уравнения do 1 i=1,nd do 1 j=1,nd ker(i,j)=gf(i,j)*dq(j) 1 continue c вычисление внеинтегрального члена call summa(nd, gf, fin, dx, f0) c переприсвоение do 2 i=1,nd fpr(i)=f0(i) 2 continue НАЧАЛО НЕЯВНОГО ЦИКЛА c IND= 100 continue call summa(nd, ker, fpr, dx, fsu) do 3 k=1,nd fou(k)=f0(k)-fsu(k) 3 continue критерий останова c po=0.


do 4 j=1,nd po=po+(fou(j)-fpr(j))** 4 continue pog=sqrt(po) решение вопроса о выходе из цикла c if(pog.lt.pogr) go to переприсвоение c do 5 j=1,nd fpr(j)=fou(j) 5 continue IND=IND+ GO TO конец неявного цикла c 6 continue write(2,40) ind write(2,40) nd dl=br-bl write(2,41) bl, br, dx, p1, p write(2,42) 42 format(5x"----------СЕТКА----------") write(2,41) (xp(i), i=1,nd) write(2,43) 43 format(5x"---ОТКЛИК СИСТЕМЫ---") write(2,41) (fou(i), i=1,nd) Stop End 4.2.4. Некоторые результаты тестовых расчтов Для возможности применения программы численного решения интеграль ного эволюционного уравнения е необходимо протестировать на примерах, имеющих аналитическое решение.

В качестве тестового примера приведм решение задачи Коши для линейно го дифференциального уравнения первого порядка вида dy 3 y 3, (4.2.13) dt t t с начальным условием ya 0.

Общее решение однородного уравнения, соответствующего уравнению (4.2.13), легко получается разделением переменных и последующим интегрирова C нием: y t. Применяя метод вариации произвольной постоянной, получаем t общее решение неоднородного уравнения (4.2.13):

2 A y t 2 3, t t где A – произвольная постоянная, значение которой определяется из начального условия: A 2a. Таким образом, решение поставленной задачи Коши имеет вид:

2 2a y t 2 3.

t t Для сведения поставленной задачи Коши к интегральному уравнению пре образуем коэффициент при y в уравнении (4.2.13) следующим образом:

3 2 2 2 2.

t tt Уравнение (4.2.13) принимает вид 3 dy 2 y 2 y 3. (4.2.14) t dt t Преобразованное уравнение (4.2.14) имеет вид (2.8), где p0 2, pt 2, f t 3 (4.2.15) t t и, следовательно, для него применим алгоритм решения, рассмотренный во вто рой главе. Нетрудно видеть, что интегральное уравнение, к которому сводится поставленная задача Коши, и функция Грина имеют, соответственно, вид 3 b t 0 t Gt, 2 d, (4.2.16) a b d, Gt, exp 2t, Gt, 0 при t a и t.

где 0 t G t, a Для выбранного тестового примера численное решение основной задачи – уравнения (4.2.16), находилось при помощи программы Echo, реализующей раз ностный алгоритм, при значениях начальной и конечной точек промежутка a, b и коэффициента для модели нулевого приближения a 1, b 6, p0 2. Функ циональные зависимости для “управляющего параметра” и входного сигнала име 2, f t 3.

ли вид (4.2.15): pt t t Для того же тестового примера численное решение основной задачи – урав нения (4.2.16), находилось для тех же значений параметров при помощи програм мы Direct, реализующей алгоритм последовательных приближений.

Собственные значения ядра интегрального уравнения, как правило, неиз вестны [25]. Поэтому для обнаружения и исключения случая “попадания” в окре стность собственного значения все расчеты проводились на последовательно из мельчаемых разбиениях. Результаты аналитического и численного решения при помощи разностного алгоритма (программа Echo) приведены ниже в таблице 4.2. только для мощности разбиения 100. Результаты численного решения посредст вом алгоритма последовательных приближений (программа Direct) в таблице не отражены. Графики аналитического решения, численного решения при помощи программы Echo и при помощи программы Direct (для значений погрешности 0, и 0,01) приведены на рисунках 4.2.3 и 4.2.4.

Таблица 4.2.1. Сравнение аналитического и численного решения при помощи разностного алгоритма основной задачи Аналитическое решение Численное решение 100 Сетка: 1.000 6.000 5.000.050 2. 1.025 1.075 1.125 1.175 1.225 ----------СЕТКА--------- 1.275 1.325 1.375 1.425 1.475 1.025 1.075 1.125 1.175 1. 1.525 1.575 1.625 1.675 1.725 1.275 1.325 1.375 1.425 1. 1.775 1.825 1.875 1.925 1.975 1.525 1.575 1.625 1.675 1. 2.025 2.075 2.125 2.175 2.225 1.775 1.825 1.875 1.925 1. 2.275 2.325 2.375 2.425 2.475 2.025 2.075 2.125 2.175 2. 2.525 2.575 2.625 2.675 2.725 2.275 2.325 2.375 2.425 2. 2.775 2.825 2.875 2.925 2.975 2.525 2.575 2.625 2.675 2. 3.025 3.075 3.125 3.175 3.225 2.775 2.825 2.875 2.925 2. 3.275 3.325 3.375 3.425 3.475 3.025 3.075 3.125 3.175 3. 3.525 3.575 3.625 3.675 3.725 3.275 3.325 3.375 3.425 3. 3.775 3.825 3.875 3.925 3.975 3.525 3.575 3.625 3.675 3. 4.025 4.075 4.125 4.175 4.225 3.775 3.825 3.875 3.925 3. 4.275 4.325 4.375 4.425 4.475 4.025 4.075 4.125 4.175 4. 4.525 4.575 4.625 4.675 4.725 4.275 4.325 4.375 4.425 4. 4.775 4.825 4.875 4.925 4.975 4.525 4.575 4.625 4.675 4. 5.025 5.075 5.125 5.175 5.225 4.775 4.825 4.875 4.925 4. 5.275 5.325 5.375 5.425 5.475 5.025 5.075 5.125 5.175 5. 5.525 5.575 5.625 5.675 5.725 5.275 5.325 5.375 5.425 5. 5.775 5.825 5.875 5.925 5.975 5.525 5.575 5.625 5.675 5. Отклик системы: 5.775 5.825 5.875 5.925 5. ----------ОТКЛИК СИСТЕМЫ--------- 0.046 0.121 0.176 0.216 0. 0.265 0.279 0.289 0.294 0.296.089.155.203.239. 0.296 0.294 0.291 0.287 0.282.283.295.303.307. 0.277 0.271 0.265 0.259 0.253.308.305.302.298. 0.247 0.241 0.234 0.228 0.222.287.281.275.269. 0.217 0.211 0.205 0.200 0.195.256.250.244.238. 0.189 0.184 0.180 0.175 0.170.226.220.215.209. 0.166 0.162 0.158 0.154 0.150.199.194.189.184. 0.146 0.143 0.139 0.136 0.133.175.171.167.163. 0.130 0.126 0.124 0.121 0.118.156.152.149.145. 0.115 0.113 0.110 0.108 0.105.139.136.133.130. 0.103 0.101 0.099 0.097 0.095.125.122.119.117. 0.093 0.091 0.089 0.087 0.086.112.110.108.106. 0.084 0.082 0.081 0.079 0.078.102.100.098.096. 0.076 0.075 0.073 0.072 0.071.093.091.090.088. 0.069 0.068 0.067 0.066 0.065.085.084.082.081. 0.063 0.062 0.061 0.060 0.059.078.077.076.075. 0.058 0.057 0.056 0.055 0.055.072.071.070.069. 0.054 0.053 0.052 0.051 0.050.067.066.065.064. 0.050 0.049 0.048 0.047 0.047.063.062.061.060..058.058.057.056. Мощность разбиения результат 0, аналитического решения 0, 0, результат 0, численного решения 0, 0, результат решения 0, итерационным методом.

0, погрешность 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, Рис. 4.2.3. Результаты численного моделирования (программа Echo) Мощность разбиения результат аналитического 0, решения 0, 0, результат 0,200 численного решения 0, 0, результат 0,050 решения итерационным 0,000 методом.

погрешность 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, Рис. 4.2.4. Результаты численного моделирования (программа Direct) Из таблицы и графиков видно, что для сравнительно небольшой мощности разбиения промежутка времени эволюции оба использованных алгоритма чис ленного решения дают результаты, близкие к аналитическому решению. Для ите рационного алгоритма точность вычислений растт с уменьшением задаваемой погрешности расчтов, что находится в полном соответствии с доказанной выше абсолютной и равномерной сходимостью борновского ряда.

4.2.5. Результаты численного решения основной задачи математического моделирования динамики малого предприятия разностным методом и методом последовательных приближений В работах [20, 67] эволюция малого предприятия, как хозяйствующей еди ницы, описана в рамках динамической модели Кейнса, то есть, уравнением вида dyt pt y(t ) I (t ). (4.2.17) at В уравнении (4.2.17) переменный коэффициент pt имеет вид 1 c 1 t f pt, (4.2.18) 1 2 k 1 t где изменн знак в числителе. Входящие в коэффициент (4.2.18) величины имеют следующий смысл (см. также введение):

1 и 2 – ставки налогообложения на объем выпуска и прибыль соответст венно;

с – удельная себестоимость выпуска продукции;

k – коэффициент, отражающий долю реинвестируемых средств прибыли, не имеющих льгот по налогообложению, и оцениваемый статистическим путем 0 k 1;

t – доля чистой прибыли, отчисляемой на реинвестирование 0 t 1 ;

f – показатель фондоотдачи;

I t – внешние инвестиции, получаемые предприятием на безвозмездной основе.

Для фоновой модели с постоянной долей чистой прибыли, отчисляемой на реинвестирование, для постоянного коэффициента p0 примем значение 1 c 1 0 f p0.

1 2 k 1 Дифференциальное эволюционное уравнение имеет вид dyt p0 y(t ) pt y(t ) I t, at где 1 c 1 A t 1 c 1 0 f.

pt pt p0 b2 A 2 1 2 k 1 1 k 1 t b Для доли чистой прибыли, отчисляемой на реинвестирование, примем функцио нальную зависимость от времени вида [28]:

t 2, t b, t A, t b.

Здесь A максимально допустимая доля чистой прибыли, отчисляемая на реинве стирование. Так как A b b 2, то для темпа реинвестирования получаем A величину.

b Результаты численного решения эквивалентного основной задаче (задаче Коши с нулевым начальным условием) интегрального эволюционного уравнения методом последовательных приближений для значений параметров a 1, b 6, 0,0001t p0 0,00009, pt 0,00009, I t 5t приведены на 1 0,0221 0,013t рисунке 4.2.5.

Модель малого предприятия 80000, 70000, 60000, 50000, 40000,000 Результат численного решения 30000, 20000, 10000, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, Рис. 4.2.5. Результаты численного моделирования динамики малого предприятия Как видно из графика, вначале рассматриваемого периода, идет резкое уве личение стоимости основных производственных фондов за счет собственных средств и внешних инвестиций. Несмотря на возрастание внешних инвестиций, тем не менее, со временем происходит постепенное уменьшение стоимости ос новных средств, связанное с их износом.

4.3. Примеры численного решения основной задачи для многомерной модели параметрической системы 4.3.1. Трхмерная модель параметрической системы В общем случае [27, 64] (или параграф 1.3) динамика трхмерной экономи ческой системы с непрерывным временем описывается системой из трх обыкно венных дифференциальных уравнений в нормальной форме dx dt p1 t x p2 t x p3 t x f, 1 1 1 2 1 3 dx p12 t x1 p2 t x 2 p3 t x 3 f 2, 2 (4.3.1) dt dx p 3 t x1 p 3 t x 2 p 3 t x 3 f 3, dt 1 2 коэффициенты которой имеют экономический смысл прямых материальных за трат и, вообще говоря, зависят от времени неявно, через параметры, аналогично тому, как зависит коэффициент (В.11) в уравнении (В.10), описывающем динами ку малого предприятия.


Установление для коэффициентов прямых материальных затрат функцио нальных зависимостей вида p ij t p ij c ij, f ji, 1, 2 (4.3.2) требует проведения длительных экспериментальных исследований, которые должны проводиться для отдельно взятых классов однотипных экономических систем (см. введение). Для того чтобы зависимости вида (4.3.2) достаточно хоро шо приближали экономическую реальность, их необходимо проводить в течение длительного времени, учитывая различные факторы, влияющие на динамику. По добные исследования проводятся методами эконометрики [28]. Применение ме тодов эконометрики к исследованию конкретных экономических систем не явля ется, однако, задачей данной работы. По указанной причине в качестве примеров численного моделирования трхмерной параметрической системы рассмотрены три в определнной степени абстрактные параметрические модели в предположе нии, что их динамика описывается системой обыкновенных дифференциальных уравнений с коэффициентами, зависящими от времени явно.

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

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

dx e 0, 01t 2 e 0, 5t 3 t x 0,1 2 x 0,2 2 x 0,4 2 e 0, 02 t ( t 6 ) t 1 t dt 20 dx 1 2 0,1t 0,1t ( t 6 ) 0, 06 t e e e t 0,2 2 x 0,1 2 x 0, (t 1) 2 x 30 2 2 (t 1) (t 1) dt dx 0, 05 t ( t 6 ) 0, 4 t t2 e 0,3 0,1 2 x x 0,1e x t ( t 1) 2 (t 1) 0,1 e 2, dt 0,6t 20 10 (4.3.3) Придавая системе (4.3.3) “экономический смысл”, можем считать, что координа ты вектора переменных имеют следующий смысл:

x1 – выручка;

х2 – чистая прибыль;

х3 – активы (все, что есть у предприятия, чем оно может распоряжаться).

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

Итак, решаем задачу Коши:

dx1 t dt 0,4 x 0,1x 0,2 x 20 10, 1 2 x 1 (0) 0, dx t x (0) 0, 0,2 x1 0,1x 2 0,2 x 3, (4.3.4) dt3 30 10 x 3 ( 0) 0.

t dx 3 0,1x 0,1x 0,1x t, 1 2 dt 20 10 Полагая x 1 1 e t, x 2 2 e t, x 3 3 e t, получаем:

dx1 dx 2 dx 1 ( )e t, 2 e t, 3 e t.

dt dt dt Подставляя производные в однородную систему уравнений, получаем СЛАУ (0,4 ) 1 0,1 2 0,2 3 0, 0,2 (0,1 ) 0,2 0, 1 2 0,1 1 0,1 2 (0,1 ) 3 0, характеристическое уравнение которой имеет вид:

0,4 0,1 0, 0,1 0,2 0, 3 0,6 2 0,11 0,006 0.

0, 0,1 0, 0, Решая характеристическое уравнение, получаем собственные значения:

1 0,1;

2 0,2;

3 0,3.

Находим собственные векторы:

для 1 0,1 собственный вектор x1 1 ;

для 2 0,2 собственный вектор x2 0 ;

для 3 0,3 собственный вектор x3 1.

Записываем фундаментальную матрицу системы решений однородной системы уравне ний e 0,1t e 0,3t e 0, 2t X (t ) e 0,1t e 0, 3t.

e 0,1t e 0, 2t Находим матрицу X 1 (t ) :

det X (t ) e 0, 6t e 0, 6t e 0, 6t e 0, 6t, e 0,5t e 0,5t e 0,1t e 0,1t e 0,5t e 0,1t 1 X 1 (t ) 0, 6t e 0, 4t e 0, 4t 0 e 0, 2t e 0, 2t 0, e 0, 3t e 0, 3t e 0, 3t e 0, 3t e 0 1 1 X 1 (0) 1 1 0.

1 0 Матричная функция Грина дифференциального оператора фоновой системы (4.3.4) в со ответствии с формулой (1.3.10) принимает вид:

e 0,1t e 0,3t 1 e 0, 2t G (t ) H (t ) X (t ) X 1 (0) H (t ) e 0,1t 1 1 0, 3t 0 e e 0,1t 0 1 0 e 0, 2 t e 0,1t e 0, 2t e 0,3t e 0,1t e 0,3t e 0,1t e 0, 2t H (t ) e 0,1t e 0,3t e e.

e 0,1t 0,1t 0, 3t (4.3.5) e 0,1t e 0, 2t e 0,1t 0, 2t e 0,1t e Решение задачи Коши по принципу Дюамеля записывается в виде интеграла z z z z e 0,1( t ) e 0, 3( t ) 20 0,1 e 0,1( t ) e 0, 2 ( t ) e 0, 3( t ) e 0,1( t ) e 0, 2 ( t ) t 3 0 0,2 d.

e 0,1( t ) e 0, 3( t ) e e 0,1( t ) 0,1( t ) 0, 3 ( t ) e 20 0,3 0, e 0,1( t ) e 0, 2 ( t ) e 0,1( t ) e 0, 2 ( t ) e 0,1( t ) Для первой координаты имеем:

d.

e e 0,3(t ) t z 1 t e 0,1( t ) 2 0, 2 ( t ) 20 19 60 0, 2 2 0 7 60 0, Найдем каждое из трех слагаемых отдельно, используя метод интегрирования по частям:

t 2 t 401 401 0,1t 2 2 0 19 6 0 0,2 d 0,1( t ) t e e, 26 3 t d t 11 11 0, 2 t 0, 2 ( t ) e e 60 0,1, 12 12 t 2 41 205 205 0,3t t 20 7 20 d 0, 3( t ) 2 t e e.

6 18 27 Окончательно получаем, что 401 0,1t 11 0, 2t 205 0,3t 1 2 389 z1 e t t e e.

3 12 27 3 36 Аналогично находим другие координаты.

d, e t z 2 t e 0,1( t ) 2 0, 3 ( t ) 2 0 19 6 0 0, 2 2 0 7 2 401 0,1t 205 0,3t 1 2 98 z2 t t e e ;

3 27 3 9 e t d, z 3 (t ) e 0,1( t ) 2 0, 2 ( t ) 20 19 60 0, 2 6 0 0, 401 0,1t 11 0, 2t 1 2 157 z3 e t t e.

3 12 2 12 Таким образом, решение задачи Коши имеет вид:

401 0,1t 11 0, 2t 205 0,3t 1 2 389 x1 e t t e e, 3 12 27 3 36 401 0,1t 205 0,3t 1 2 98 x2 t t e e, (4.3.6) 3 27 3 9 401 0,1t 11 0, 2t 1 2 157 x3 e t t e.

3 12 2 12 Результаты численного моделирования динамического поведения парамет рической системы (4.3.3) представлены в графическом виде на рисунках 4.3.1 – 4.3.3. Там же представлены графики результатов аналитического решения (4.3.6) задачи Коши для фоновой системы (4.3.4).

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

dx e 0, 01t 2 e 0, 5t 3 t 0,4 e 0, 02t ( t 6 ) x1 0,1 x 0, t 1 x 20 10, t 1 dt dx 2 e 1 e 2 e 0,1t 0,1t ( t 12 ) 0, 06 t t 0,2 (t 1) 2 x 0,1 (t 1) 2 x 0,1 (t 1) 2 x 30 10, (4.3.7) dt dx 0, 05 t ( t 6 ) e 0, 4t t 0,3 3 0, (t 1) 2 x 0,1 e 0, 6t 2,9 x 0,1e x t.

( t 1) 2 dt 20 10 Очевидно, что модельной системе уравнений (4.3.7) соответствует та же фоновая система (4.3.4) и, следовательно, те же функция Грина (4.3.5) и аналитическое решение (4.3.6) задачи Коши (4.3.4). Функциональные зависимости коэффициен тов системы уравнений (4.3.7) естественно существенным образом отличаются от соответствующих функциональных зависимостей системы уравнений (4.3.3).

Сравнение результатов численного моделирования для параметрической системы, описываемой системой уравнений (4.3.7), с результатами аналитическо го решения для фоновой системы (4.3.6) приведено на рисунках 4.3.4 – 4.3.6.

Модель 3. Пусть модель динамики трхмерной параметрической системы описывается системой трх обыкновенных дифференциальных уравнений (4.3.1) dy dt p1 t y p 2 t y p3 t y f, 1 1 1 2 1 3 dy p12 t y 1 p 2 t y 2 p3 t y 3 f 2, 2 (4.3.8) dt dy p 3 t y 1 p 3 t y 2 p 3 t y 3 f 3, dt 1 2 с коэффициентами 2t t2 t p1 t 0,4 0,008 e 0, 6t sin, 12 2t t 3 t 2 0, 5t p2 t 0,1 0,0009 e, cos 8 2t t 4 t 3 0, 7 t p t 0,2 0,007 e, sin 48 2t t2 t p12 t 0,2 0,05 e 0,8t cos, 8 2t t3 t p2 t 0,1 0,003 e 0, 4t sin, (4.3.9) 28 2t t 3 t 2 0, 6 t p3 t 0,2 0,004 e, cos 27 2t t2 t p t 0,1 0,001 e 0,9t sin, 18 2t t3 t p2 t 0,1 0,002 e 0,5t cos, 6 2t t4 t p3 t 0,1 0,003 e 0, 4t sin 42 и правыми частями t 0,0112t t 2, f 3 t 0,02t 2 12t 13.

f 1 t t 1, f (4.3.10) Динамика фоновой системы описывается системой обыкновенных дифференци альных уравнений (4.3.4) с той же постоянной матрицей коэффициентов. Графи ческие зависимости численного решения задачи Коши для параметрической сис темы в сравнении с графиками аналитического решения для фоновой системы, приведены на рисунках 4.3.7 – 4.3.9.

Для моделей 1 и 2 результаты численного решения задачи Коши для пара метрической системы отличаются от результатов аналитического решения задачи Коши для фоновой системы не слишком сильно в силу того, что функциональные зависимости коэффициентов соответствующих систем дифференциальных урав нений являются “слабыми”. Для модели 3, эти зависимости являются более силь ными и одновременно более сложными. Соответствующие результаты численного моделирования, сохраняя качественный вид аналитического решения для фоно вой системы, отличаются от результатов аналитического решения в значительно большей степени.

X 0, 0, 0, 0, 0,08 Числ.метод Аналитика 0, 0, 0, 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1 1,2 t Ри с. 4.3.1. Численное моделирование динамического поведения параметрической модели 1, X X 0, 0, 0, Аналитика Числ.метод 0, 0, 0 0,1 0,2 0,3 0,4 0,5 0,60,7 0,8 0,9 1,0 1,1 1,2 t Рис. 4.3.2. Численное моделирование динамического поведения параметрической модели 1, X X 0, 0, 0, 0 0,1 0,2 0,30,4 0,5 0,6 0,7 0,80,9 1,0 1,1 1,2 t -0, -0, Аналитика Числ.метод -0, -0, -0, -0, -0, Ри с. 4.3.3. Численное моделирование динамического поведения параметрической модели 1, X X 0, 0, 0, 0, 0,08 Числ.метод Аналитика 0, 0, 0, 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1 1,2 t Рис. 4.3.4. Численное моделирование динамического поведения параметрической модели 2, X X 0, 0, 0, Аналитика Числ.метод 0, 0, 0 0,10,20,30,40,50,60,70,80,91,01,11,2 t Рис. 4.3.5. Численное моделирование динамического поведения параметрической модели 2, X X 0, 0, 0, 0 0,10,20,30,40,50,60,70,80,91,01,11,2 t -0, -0, Аналитика Числ.метод -0, -0, -0, -0, -0, Рис. 4.3.6. Численное моделирование динамического поведения параметрической модели 2, X X 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 t -0, - Числ.метод -1, Аналитика - -2, - Рис. 4.3.7. Численное моделирование динамического поведения параметрической модели 3, X X 0, 0, 0, 0, Аналитика Числ.метод 0, 0, 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 t Рис. 4.3.8. Численное моделирование динамического поведения параметрической модели 3, X X 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 t -0, -0, -0, -0, Аналитика Числ.метод -0, -0, -0, -0, Рис.

4.3.9. Численное моделирование динамического поведения параметрической модели 3, X 4.4. Пример численного моделирования переходного процесса в параметрической модели ИИП 4.4.1. Структурные параметры модели ИИП В параграфе приведн пример численного моделирования динамики индук тивного измерительного преобразователя в рамках предложенной параметриче ской математической модели (параграфы 1.2, 2.3).

Положим, что световой день, в течение которого проводятся измерения, длится 10 часов: с 7-00 часов утра до 17-00 часов вечера. Элементы a ij фоновой матрицы A определяются е выражением (1.2.11) из параграфа 1.2:

r A L L.

1 C RC Здесь значения сопротивления, индуктивности и емкости цепи выбраны в отсчт ный момент времени наблюдений – 7 часов утра (для удобства вводится в про грамму как t0 0 ). Вычисление изменения параметров за время релаксации пере ходного процесса осуществляются по формулам (1.2.17) – (1.2.17):

1 dRT RT RT0 1 T T T T0, RT0 dT 1 dLT LT LT0 1 T T T T0, LT0 dT 1 dCT C T C T0 1 T T0, C T0 dT T T где R0, L0, C0, r0 – фоновые значения параметров цепи в начале светового дня, то есть в момент t0 0. Значения R0, L0, C0, r0 для идеальной модели взяты из [10]:

R0 4,5 кОм, L0 67 мГн, C0 3,9 нФ, r0 175 Ом Электродвижущая сила, вызванная в индуктивном измерительном преобразовате ле изменением первичного поля, вычисляется по формуле [10] Et E0 t 1t t ф, где E0 300 мкВ – значение ЭДС в начальный момент времени, 1t, 1 t t ф – единичные скачки, t ф 104 с – время релаксации переходного процесса.

Для упрощения расчтов зависимость температуры от времени выберем в виде функции “шапочка” [12] (рисунок 4.4.1), достаточно реалистично описы вающей ход температуры в течение светового дня в средней полосе России:

2, a t a, B exp t a T t 0, t a, t a.

Здесь B 923.639, 71196, a 18000.

Рис. 4.4.1. Зависимость температуры от времени в течение светового дня 4.4.2. Результаты численного моделирования Для численного моделирования введм разбиение оси времени (рисунок 4.4.2). Диапазон изменения времени t 0, 36000 с., начальный момент времени t0 0 соответствует 7 часам утра, конечный момент времени t 36000 с. соот ветствует 17-00 ч вечера.

Для упрощения расчтов предположим, что температура в 7 часов утра и в 17 часов вечера принимает равные значения 10 С, максимум температуры дости гается в полдень 12 часов и составляет 30 С. Такое предположение является чис то “техническим” и не является ограничительным. При расчетах температура пе реводится в систему СИ (градусы Кельвина).

Время переходного процесса t, c 10- t, c 0 t, ч 7 12 Рис. 4.4.2. Разбиение оси времени в течение светового дня Предполагается, что измерения проводятся каждый час. В результате чис ленных расчтов в каждом из отрезков t 10, отсчитываемых от начальных моментов времени t0, t1, t2,, t10, получены значения силы тока и напряжения.

Результаты расчтов выборочно для занеснных в таблицу 4.4.1 значений времени представлены на рисунках 4.4.3 – 4.4.6.

Таблица 4.4.1. Начальные значения для частичных промежутков времени релаксации переходного процесса Точки t0 t3 t6 t разбиения Время в 7 10 13 часах Время в 0 10800 21600 секундах Толстая линия – численный расчт для параметрической модели, тонкая линия – аналитический расчт для идеальной модели.

Время в часах t 7. Время в секундах t 0.

Рис. 4.4.3. Результаты численного моделирования для начального момента t Время в часах t 10. Время в секундах t 10800.

Рис. 4.4.4. Результаты численного моделирования для начального момента t Время в часах t 13. Время в секундах t 21600.

Рис. 4.4.5. Результаты численного моделирования для начального момента t Время в часах t 16. Время в секундах t 32400.

Рис. 4.4.6. Результаты численного моделирования для начального момента t Из представленных графиков видно, что динамика параметрической и иде альной систем существенно различается.

ЗАКЛЮЧЕНИЕ На основе проведнного в диссертационной работе анализа распро страннности параметрических систем в практических исследованиях, сделан вывод о том что, структурные параметры подавляющего числа предметных систем, взаимодействующих с окружающей средой, вследствие внешних воздействий (возмущений), возможно, опосредовано через внешние пере менные, например, температуру, зависят от времени и, следовательно, систе мы являются параметрическими.

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

Создан комплекс вычислительных программ на языке ФОРТРАН, эф фективность которого продемонстрирована на модельных примерах. При помощи разработанных вычислительных программ прямым численным мо делированием показано, что динамика линейных параметрических систем существенно отличается от динамики идеальных систем. Это отличие тем существеннее, чем сильнее внешняя среда воздействует на систему. Такое поведение систем, взаимодействующих с внешней средой, позволяет реко мендовать рассматривать их в процессе различных исследований как пара метрические системы. Эти рекомендации существенны, по крайней мере, для рассмотренных в работе предметных областей, а именно, для полевых изме рительных систем в геофизике и для динамических систем с непрерывным временем в экономике.

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

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

СПИСОК ЛИТЕРАТУРЫ 1. Адомиан Дж. Стохастические системы/Дж. Адомиан. – М.: Мир, 1987. – 376 с.

2. Астраханцев Г. В. Индукционное зондирование при изучении контрастных по электропроводности сред/Г. В. Астраханцев. – Свердловск: Институт гео физики УрО АН СССР, 1988. – 182 с.

3. Афонский А. А. Измерительные приборы и массовые электронные измере ния/ А. А. Афонский, В. П. Дьяконов. – М.: СОЛОН-ПРЕСС, 2007. – 544 с.

4. Бакалов В. П. Основы анализа цепей/В. П. Бакалов, О. Б. Журавлва, Б. И.

Крук. – М.: Горячая линия – Телеком, 2007. – 591 с.

5. Бардзокас Д. И. Математическое моделирование в задачах механики свя занных полей. Введение в теорию термопьезоэлектрическтва. Том 1/Д. И.

Бардзокас, А. И. Зобнин, Н. А. Сенник, М. Л. Фильштинский. – М.: ЕДИТО РЕАЛ-УРРС, 2005. – 312 с.

6. Бендат Дж. Прикладной анализ случайных данных/ Дж. Бендат, А. Пирсол.

– М.: Мир, 1989. – 540 с.

7. Бессонов Л. А. Теоретические основы электротехники. Электрические це пи/Л. А. Бессонов. – М.: ГАРДАРИКИ, 2007. – 7091 с.

8. Брке У. Пространство-время, геометрия, космология/ У Брке. – М.: Мир, 1985. – 416 с.

9. Богословский В. А. Геофизика/В. А. Богословский, Ю. И. Горбачв, А. Д.

Жигалин, В. К. Хмелевской. – М.: КДУ, 2008. – 320 с.

10. Вакульский А. А. Аппаратура для аэрогеофизической разведки с магнит ным и электромагнитным информационными каналами/ А. А. Вакульский, Л.

Я. Мизюк и др. – Киев: НАУКОВА ДУМКА, 1985. – 256 с.

11. Васильев Д. В. Радиотехнические цепи и сигналы/Д. В. Васильев, М. Р.

Витоль, Ю. Н. Горщенков и др. – М.: РАДИО И СВЯЗЬ, 1982. – 528 с.

12. Владимиров В. С. Уравнения математической физики/В. С. Владимиров.

– М.: Наука, 1981. – 512 с.

13. Волкова В. Н. Теория систем/В. Н. Волкова, А. А. Денисов. – М.: Высшая школа, 2006. – 511 с.

14. Гантмахер Ф. Р. Теория матриц/Ф. Р. Гантмахер. – М.: Наука, 1988. – с.



Pages:     | 1 | 2 || 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.