Алгоритм численного решения задачи Шоуолтера-Сидорова для систем леонтьевского типа Текст научной статьи по специальности «Математика»

УДК 519.622.2

АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ ЗАДАЧИ ШОУОЛТЕРА-СИДОРОВА ДЛЯ СИСТЕМ ЛЕОНТЬЕВСКОГО ТИПА

А.В. Келлер

ALGORITHM OF NUMERICAL SOLUTION OF SHOULTER-SIDOROV PROBLEM FOR LEONTIEV-TYPE SYSTEMS

A.V. Keller

В статье рассмотрено решение задачи Шоуолтера-Сидорова для систем леонтьев-ского типа с необратимым оператором при производной. Рассмотрение начальных данных Шоуолтера-Сидорова позволяет расширить спектр практического применения модели. Предложен алгоритм ее численного решения, исследована сходимость.

Ключевые слова; задача Шоуолтера-Сидорова, системы леонтьевского типа, исследование сходимости.

The article deals with the solution of Shoulter-Sidorov problem for Leontiev-type systems with irreversible operator at the derivative. Consideration of the initial data of Shoulter-Sidorov enables greater range of practical applications of the model. An algorithm for its numerical solution has been proposed, the convergence has been investigated.

Keywords: Shoulter-Sidorov problem, Leontiev-type systems, convergence investigation.

Введение

В [1, 2] предложен численный алгоритм решения задачи Коши

и(0) = и0 (1)

для линейной системы обыкновенных дифференциальных уравнений

й = Su + g, (2)

основанный на идеях теории полугрупп операторов. Здесь 5 - квадратная матрица порядка п. В [3, 4] этот подход был распространен на задачу (1) для вырожденной системы уравнений

Ы = Ми + / (3)

с использованием идей теории вырожденных полугрупп операторов [5] (здесь Ь и М - квадратные матрицы, порядка п, причем Ь = 0). Одним из важных случаев системы (3) является хорошо известная система В.В. Леонтьева «затраты-выпуск» с учетом запасов (см. в [6]), поэтому в [3] было предложено такие системы уравнений называть «системами леонтьевского типа».

Простота предложенного в [3, 4] алгоритма обеспечивает высокое качество получаемого программного продукта, что выгодно отличает данный алгоритм от использовавшихся ранее методов Эйлера, Рунге-Кутта, итерационных и других методов [7-9]). Основным недостатком этого алгоритма (как впрочем, и всех остальных)

является принципиальная неразрешимость задачи (1) для системы (3) при произвольных начальных векторах щ е . Эта трудность преодолевается, например в [3], где вектор-функция / :[0,Т]-»91" постоянная, введением множества

допустимых начальных значений, понимаемых как фазовое пространство системы (3). В [4] при снятии условия постоянства вектор-функции / налагаются условия согласования / с начальным значением м0. Заметим, что условия согласования в том или ином виде имеют место и во всех других алгоритмах.

Для условий согласования (как и для построения фазового пространства) необходимы проекторы, которые либо выражаются через контурные интегралы от матриц-функций, либо являются пределами матричных последовательностей. Ввиду неустойчивости любого проектора относительно малых возмущений такое вычисление матрицы проектора очень затруднительно. Поэтому, например в [10], при построении системы (3), моделирующей экономику коммунального хозяйства, пришлось ограничиться малыми городами, т.е. такими, где матрицы Ь и М имеют порядок не больше 10. Именно малость порядка матриц Ь и М сделало возможным вычисление проекторов «вручную».

Келлер Алевтина Викторовна - канд. физ.-мат. наук, Keller Alevtina Viktorovna - PhD, associate professor,

доцент, зав. кафедрой общеобразовательных дисциплин head of general disciplines department of SUSU;

ЮУрГУ; [email protected] [email protected]

Между тем в современной математической литературе существуют попытки теоретического осмысления так называемых «неклассических» задач для системы (3) [11, 12], основным достоинством некоторых является однозначная разрешимость при любых начальных данных щ е 91” . Разработка численных алгоритмов решения таких задач позволит избавиться как от трудоемкого построения фазового пространства (и не менее трудоемкой редукции системы (3) к системе (2), заданной на нем), так и от трудоемкой проверки условий согласования. Основная цель данной статьи - построение алгоритма численного решения задачи Шоуолтера-Сидорова

(аЬ-М) 1Ь (и(0)-и0) = 0

для системы (3) (числа а и р вычисляются на первом шаге алгоритма). Статья кроме введения и списка литературы содержит две части. В первой дается теоретическое обоснование алгоритма, а во второй приведены основные этапы вычислений.

1. Задача Шоуолтера-Сидорова

Пусть Ь и М - квадратные матрицы порядка п, причем, следуя [13], гл. XII, п. 2, пучок матриц /лЬ-М назовем регулярным, если существует число Я е С такое, что М} Ф О .

Заметим, что условие регулярности пучка матриц эквивалентно условию Ь -регулярности матрицы М [3, 4]. Поэтому, как показано в [5], гл. 4, при условии регулярности пучка существуют единственным образом определяемые матрицы Я, 5,

М0, X], Q порядка п, такие, что Ь-резольвента (рЬ-М)~1 матрицы М разлагается в ряд Лорана

(рЬ-М)-1^ //^0(7-6) +

1—0

(4)

+2>_/ям4е

1=1

в окрестности бесконечно удаленной точки, причем Н - нильпотентная матрица со степенью нильпотентности р, Q - идемпотентная матрица, ММ0, М0М, ЦЬ и Ых - диагональные матрицы с нулями и единицами на главной диагонали. Поскольку с1е1:(АХ-М)^0, то многочлен с!е1:(/1/,-А/) = 0 имеет не более п различных нулей, которые расположены в круге радиуса а, а значит, при \/и\ > а разложение (4) имеет

место. Точка оо называется устранимой особой точкой Ь -резольвенты матрицы М, если р = 0 в (4); и полюсом порядка р е N в противном случае. В дальнейшем, немного отходя от клас-

сического стандарта, будем называть устранимую особую точку полюсом порядка нуль. Итак, пусть пучок /лЬ-М регулярен, И 00 - полюс

порядка р е {0}иК; тогда можно выбрать число

а и рассмотреть задачу Шоуолтера-Сидорова

[^(М)]/,(М(0)-Мо) = 0, (5)

для системы уравнений леонтьевского типа

Хм = Ми + /, (6)

где К1а(М) = [аЬ-М) 1 Ь -правая Ь-резольвента матрицы М, в отличие от ее левой Ь -резольвенты 11а(М) = Ь{аЬ-Му1, а /:[О,Т] -»Шп -некоторая вектор-функция.

Решением системы (3) называется вектор-функция щ еС1((0;Т);ЭТ")пс([0;Т];91"),

удовлетворяющая уравнениям системы. Решение системы (6) называется решением задачи (5), (6), если оно вдобавок удовлетворяет уравнениям (5). Имеет место [5, гл. 4]

ТЕОРЕМА 1.1. Пусть пучок /лЬ-М регулярен, /?е{0}кЖ - порядок полюса Ь -резольвенты матрицы М, вектор-функция f: [0,Т] -> 9?"

такая, что {ъьа(М)^ /е с([0;Т];9?” |, а

1-(4(АО)Р/бС^+1((0;Т);9Г)пС'([0;Т];Я").

Тогда при любом и0 е 9?” существует единственное решение задачи (5), (6), которое к тому же имеет вид

к(0 = -Е нЧм01 (I - <2)/(9)(0+и\ +

д=0

t

о

Здесь

и'

г

* ‘УгА^-мТ

Г

г

контур у = {// 6 С: 1^1 = г > а|. Контурные интегралы не очень удобны в численных расчетах, поэтому в [3, 4] предложен другой подход, основанный на аппроксимациях типа Уиддера-Поста [5, гл. 2]. Именно справедлива

ТЕОРЕМА 1.2. Пусть пучок цЬ-М регулярен, - порядок полюса Ь -резоль-

венты матрицы М. Тогда

1ІШ

к->°°

Vі

*(/7+1)

L —

к(р +1)

-М

= и‘,

limj_kLLk(M)f+1) = Q,

k->aо

lim

к-> оо

і—

*Ср+і)

-м

х—

-м

R‘.

к(р +1)

Теперь пусть пучок /лЬ-М регулярен, р £ | О | ^ N - порядок полюса X -резольвенты матрицы М в точке оо. Фиксируем Т е Ш+, ?е(0,Т), £eN и положим

л-i

к(р+1)

X —

Кр+1)

п(р+1)

-М

Qk=[kLLk(M)]P ,

I X-------:—M

1 КР +1)

ч-l

-M

k(p+1)

Выберем вектор u0 є SR”, вектор-функцию / є Cp+1 ((0; T); 5R") n Cp ([О; T]; 9Г ) и построим вектор-функцию

«* (0=-І н9мо' С1 - а)/(?) (о+икЧ+

?=О

+/#-'&/(*)*■

о

ТЕОРЕМА 1.3. Пусть пучок /лЬ-М регулярен, р є {О} - порядок полюса X-

резольвенты матрицы М в точке оо. Тогда существует константа С = С(Ь, М, Т) є 5Й+

(2

такая, что ||м(/) — (?)[| < — при всех І є [0;Т],

щ є 9Г, и f є Ср+1 ((0; Т); Я”) г» Ср ([0; Т]; Ж").

Доказательство теоремы основывается на оценках

[^(М)](р+1)-фХ

р+1

КС

р+1

/Й (.Р + \)/Лк~Х РР+Х~к Р

Ы-иЛ^1^2

•• к 11 2 рр-1к взятых из [5, гл. 2], где /? є 9?+

2. Алгоритм решения

Построение алгоритма начнем с допущения, что det М 9*0. Это допущение не ограничивает общности предыдущих рассуждений. Действительно, при условии регулярности пучка /лЬ-М, можно сделать замену и = eXtv в уравнении (3) и перейти к уравнению

Lv = (М - AX)v+e~Xt f (7)

того же вида, что и (3), но det(M - XL) ф 0. Обратный переход от решений системы (7) к решениям системы (3) очевиден.

На первом шаге алгоритма нужно найти числа а е3? и pe(0}uN. Можно разумеется,

разложить X -резольвенту матрицы М в ряд (4) и тем самым сразу же найти эти числа. Однако, существует другой менее трудоемкий путь. Рассмотрим многочлен

det(yuZ - М) = ап/лп + ап_хцп~х + . + ахц + а0 . Поскольку а„ = detX, то а„ = 0. Далее, коэффициент а/ есть сумма слагаемых, каждое из которых есть произведение одного из миноров порядка 1 матрицы X на число, / = 1,..., п-\, а0 = det(-M). Поэтому степень многочлена det(fiL-M) не выше rankX, т.е. ранга матрицы X. Итак,

det{/лЬ-М) = aq/iq +ач_фч~х +... + а1(и + а0, где q = deg det{p.L - M) < rank X, поэтому если взять число ae5R таким, что с л W

\а\ >тах-Ц,\а„

ы

ч/=о

то det(aX - М) ф 0, и значит, существует матрица (аЬ-Му1. Далее, считая, что матрица М обратима, представим

det(1^^Z - М) = det М det(/^M_1X -1).

Зная, что порядок полюса в точке оо резольвенты

(///-ЛГ1!) равен нулю, легко найти, что порядок полюса X -резольвенты матрицы М в точке оо равен п-ц. Итак, числа а и /? = и-д найдены.

Тогда находя значение к, с которого можно начинать считать приближенные проекторы, получим, что при

*>0 Е Ы+1

\Я I

мы не сможем оказаться даже вблизи точки X -спектра оператора М.

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

det(//(^? + \)Ь - 1М) = а^4 [лч (р +1)9 +

+а 1?9+1^9-1(/7 + 1)9-1 + — + а^п~1 р. + ^ а0,

где aq *0,q < rank і . Тогда, учитывая p = n-q при

1

k>

,| (n-q)

/=<7+1

если И < 1

K M ^n~q>

(n-q)

l=q+1

если Ы > 1

мы не сможем оказаться даже вблизи точки Ь -спектра оператора М.

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

3. Пример Леонтьева

Операторы Ь и М зададим матрицами

|% Хо 21/ > /200

х = Хоо 103/ /200 8/ /25

0 0 0

М =

/25

-V

v /15

Я

10304189

11996000

-2/

-11/

/20

-70836357

119960000

13/

/15

Если переобозначить Ь = В, М = 1-А, то матрицы В и А почти совпадут с матрицами из классического примера [6]. «Почти» означает, что элементы т22 и т23 подобраны специально с целью упростить вычисления и отличаются от приведенных в примере чисел 22/25 и -3/5 на величины

22 -252291

т22 =-

т23 =-

25

-3

11996000

1139643

(8)

5 119960000

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

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

Перейдя к расчету системы (3), найдем L -спектр оператора М crL(M) = {0,2; 2,7}. Именно для того, чтобы точки і-спектра оператора М были рациональными, сделаны поправки (8). Точки L -спектр оператора М в исходном примере иррациональны и отличаются от найденных не большее, чем на одну сотую.

Далее по формулам, приведенным в первой части статьи, построим точное и приближенное решение. Приведем точное решение и результаты счета по алгоритму без комментариев, взяв при этом в качестве / = (2/, 21, 2?) (табл. 1,2).

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

Литература

1. Павлов, Б.В. Об одном методе численного интегрирования систем обыкновенных дифференциальных уравнений / Б.В. Павлов, А.Я. По-взнер//ЖВМиМФ. - 1973. - Т. 13, № 4.- С. 1056-1059.

2. Павлов, Б.В. Численное решение систем линейных обыкновенных дифференциальных уравнений с постоянными коэффициентами / Б.В. Павлов, О.Е. Родионова //ЖВМиМФ. —1994. -Т. 34, № 4. - С. 622-627.

3. Свиридюк, Г.А. Численное решение систем уравнений леонтьевского типа /Г.А. Свиридюк, С.В. Брычев // Изв. вузов. Математика. -2003.-№8.-С. 46-52.

4. Свиридюк, Г.А. Алгоритм решения задачи Коши для вырожденных линейных систем обыкновенных дифференциальных уравнений с постоянными коэффициентами / Г.А. Свиридюк, КВ. Бурлачко //ЖВМиМФ. - 2003. -Т. 43, №11.-С. 1677-1683.

5. Sviridyuk, G.A. Linear Sobolev Type Equations and Degenerate Semi-groups of Operators / G.A. Sviridyuk, V.E. Fedorov. — Utrecht-Boston-Koln-Tokyo: VSP, 2003.

6. Леонтьев, В.В. Межотраслевая экономика/ В.В. Леонтьев. - М.: Экономика, 1997.

Таблица 1

Точное решение задачи Шоуолтера-Сидорова

t *1 *2 *3

0 1 1 0,7692307692

1/12 1,079594653 1,115125754 0,6545486955

1/6 1,211221782 1,260474368 0,5698255794

1/4 1,399261210 1,434354380 0,5156284632

1/3 1,649051640 1,634491393 0,4925503470

5/12 1,967129067 1,857877706 0,5012141312

1/2 2,361526596 2,100584518 0,5422779145

7/12 2,842149925 2,357527831 0,6164435987

2/3 3,421246553 2,622172543 0,724463383

3/4 4,113994781 2,886165457 0,867151766

5/6 4,939240310 3,138871569 1,045399350

11/12 5,920418738 3,366797082 1,260189034

1 7,086710966 3,552864395 1,512617818

Таблица 2

Приближенное решение задачи Шоуолтера-Сидорова по алгоритму

t *1 х2 *3

0 1 1 0,7692307692

1/12 1,0795957679 1,1151265428 0,6545494099

1/6 1,2112238218 1,260475827 0,5698269423

1/4 1,3992644345 1,4343566045 0,5156305665

1/3 1,6490560881 1,6344943024 0,4925530898

5/12 1,9671347889 1,8578811956 0,5012174099

1/2 2,36153371 2,1005887151 0,5422821007

7/12 2,8421584126 2,3575325354 0,6164483714

2/3 3,4212567008 2,622177848 0,7244687609

3/4 4,1140069697 2,8861709037 0,8671579463

5/6 4,9392549984 3,1388770098 1,0454063117

11/12 5,9204365728 3,3668021207 1,2601968472

1 7,0867328788 3,552868403 1,5126263028

7. Бояринцев, Ю.Е. Линейные и нелинейные алгебро-дифференциальные системы / Ю.Е. Бояринцев. - Новосибирск: Наука, 2000.

8. Чистяков, В.Ф. Избранные главы теории алгебро-дифференциальных систем / В.Ф. Чистяков, А.А. Щеглова. -Новосибирск: Наука, 2003.

9. Бояринцев, Ю.Е. Пучки матриц и алгебродифференциальные системы / Ю.Е. Бояринцев, И.В. Орлова. - Новосибирск: Наука, 2006.

10. Брычев, С.В. Исследование математической модели экономики коммунального хозяйства малых городов: дис. ... канд. физ.-мат. наук/

С.В. Брычев. - Челябинск: Челябинский гос. ун-т, 2002.

11. Свиридюк, Г.А. Задача Веригина для линейных уравнений Соболевского типа с относительно р-секториальными операторами / Г.А. Свиридюк, С.А. Загребина // Дифференц. уравнения. — 2002. - Т. 38, № 2. - С. 1646-1652.

12. Загребина, С.А. О задаче Шоуолтера-Сидорова / С.А. Загребина // Изв. вузов. Серия «Математика». - 2007. —№ 3. - С. 22-28.

13. Гантмахер, Ф.Р. Теория матриц / Ф.Р. Гантмахер. — 4-е изд. —М.: Наука, 1988. — 552 с.

Поступила в редакцию 2 октября 2007г.