ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 1, с. 83-95
УДК 519.624
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ1-*
© 2007 г. Е. Б. Кузнецов, В. Н. Микрюков
(125993 Москва, Волоколамское ш, 4, МАИ (ГТУ)) e-mail: kuznetsov@mai.ru; v-mikryukov@mail.ru Поступила в редакцию 22.06.2006 г.
На основе метода продолжения решения по параметру исследуется численное решение начальной задачи для системы дифференциально-алгебраических уравнений с запаздывающим аргументом. Получены необходимые и достаточные условия преобразования этой задачи к наилучшему аргументу, обеспечивающему наилучшую обусловленность соответствующей системы уравнений продолжения. Таким аргументом является длина дуги, отсчитываемая вдоль интегральной кривой задачи. Разработаны алгоритмы и программы численного интегрирования задачи, основанные на методах непрерывного и дискретного продолжения. Тестовые примеры демонстрируют эффективность предложенного преобразования. Библ. 17. Фиг. 4.
Ключевые слова: система дифференциально-алгебраических уравнений, запаздывающий аргумент, наилучший параметр, непрерывное и дискретное продолжения.
1. ВВЕДЕНИЕ
Рассмотрим решение начальной задачи для системы дифференциально-алгебраических уравнений с запаздывающим аргументом. Проблема заключается в том, что требуется найти решение системы уравнений
= /(У, Ух, Ут, х, хт. *), Ш (1Л)
О (у, х, *) = 0, * е[ *0, * * ], удовлетворяющее следующим начальным условиям:
У(0 = Ф(0, У(0 = ¥(0, х(о = 0, * е [*0-т, *о), У (*о) = Ф( *о) = Уо, х (*о) = у( *о) = Хо.
Здесь ф(*), у(*), у(0 - заданные непрерывные функцииУт(*) = у(* - т), Ут = У (* - т), хт(*) = х(* - т), т > 0,
/ : К3" + 2т + К", О : К" + т + К", * е К, у : К —-К", х : К —► Кт.
Точка означает дифференцирование по переменной *. Начальные условия (1.2) должны быть согласованными, т.е. в начальной точке *0 должно выполняться равенство
О (у (*о), х (*о), *о) = о.
При отсутствии запаздывания (т = 0) задача (1.1), (1.2) является задачей Коши для системы дифференциально-алгебраических уравнений, численное решение которой исследуется во многих работах (см., например, [1]). Преобразование этой задачи к наилучшему (в смысле метода продолжения решения по параметру) аргументу приводится в [2]-[4].
При отсутствии недифференциальных соотношений (вектор-функции О) задача (1.1), (1.2) является начальной задачей для системы обыкновенных дифференциальных уравнений с запазды-
Работа выполнена при финансовой поддержке РФФИ (коды проектов 06-01-00239, 06-08-00371).
вающим аргументом, численное решение которой исследуется, например, в [5]. Преобразование этой задачи к наилучшему аргументу получено в [4], [6].
В то же время численное решение начальной задачи (1.1), (1.2), по мнению авторов, не изучено.
2. НАИЛУЧШИЙ АРГУМЕНТ ЗАДАЧИ Пусть решение задачи (1.1), (1.2) задается соотношением
Г(у, х, X) = 0, F : К" + т К" + т, (2.1)
гГГ,п + т + 1 , ,
которое в евклидовом пространстве К {у, х, X} определяет единственную гладкую интегральную кривую с началом в точке
У = Уо, х = Хо, X = Хо. (2.2)
Соотношения (2.1) можно рассматривать как систему нелинейных алгебраических или трансцендентных уравнений, содержащих параметр X, который, очевидно, является аргументом задачи (1.1), (1.2).
Используя метод продолжения решения по параметру, строим кривую множества решений системы уравнений (2.1) с началом в точке (2.2). Эта кривая получается, если решать систему
(2.1) при различных значениях параметра задачи X е [Х0, X*]. Полагая, что переменные ух (X = 1, п),
х, (1 = 1, т) являются функциями параметра задачи X, уравнения продолжения строим дифференцированием соотношений (2.1) по X:
^ + ^$ + ^=0. (2.3)
дydX дхdX дх
Для построения кривой множества решений, описываемой соотношениями (2.1), систему линейных уравнений (2.3) следует разрешить относительно производных dy/dх, dх/dX, т.е. привести ее к нормальной форме Коши. При этом, учитывая начальные условия (2.2), получаем начальную задачу для системы дифференциальных уравнений
dх/dх
агат
у х
1 дг
дх.
При разрешении системы уравнений (2.3), линейных относительно производных, могут возникнуть значительные вычислительные трудности, если матрица системы является плохо обусловленной, поэтому параметр X, по которому производится дифференцирование, не всегда может оказаться удачным. Следовательно, можно поставить вопрос о выборе наилучшего в некотором смысле параметра продолжения решения системы (2.1), а значит, наилучшего аргумента задачи (1.1), (1.2).
Пусть величины у, х, X являются функциями некоторого аргумента отсчитываемого от точки (2.2), приращение которого в каждой точке интегральной кривой задачи (1.1), (1.2) представим в виде
п т
Д| = £ а АУх + £ а+п Д х, + а« + т + 1Ах. (2.4)
X = 1 1 = 1
т,« + т + 1
Здесь ак (к = 1, п + т + 1) - компоненты вектора а = (ах, ..., а« + т + х)т е К , задающего направление оси, вдоль которой отсчитывается аргумент
Для того чтобы все направления были равноправными, вектор а должен иметь постоянную длину. Пусть это будет единичный вектор, т.е.
п + т +1
а2 = £ а2 = 1. (2.5)
к = 1
Правую часть равенства (2.4) можно рассматривать как скалярное произведение вектора а на вектор приращения функций. Придавая компонентам вектора а различные значения, можно рас-
смотреть все возможные параметры продолжения, а значит, и аргументы. Так, при а = (0, ..., 0, 1)т в качестве параметра продолжения принимается переменная X.
Уравнения продолжения решения не могут быть получены дифференцированием (2.1) по параметру ц, так как эти соотношения до решения задачи (1.1), (1.2) неизвестны.
Но уравнения продолжения могут быть получены другим способом. Пусть задача (1.1), (1.2) имеет единственную гладкую интегральную кривую. Тогда системы уравнений (1.1) и (2.1) с учетом начальных условий задают одну и ту же кривую и, следовательно, в качестве уравнений продолжения решения может быть использована система (1.1). Уравнения продолжения получим, если равенство (2.4) разделим на Ац —► 0, а систему (1.1), продифференцировав функцию О по ц, перепишем с учетом равенства
. _ с/у _ /У / ^ ц
У сИ аШц' Тогда приходим к следующей системе уравнений продолжения:
п т
агуг, ц + ^^ а 1 + пХ], ц + ап + т +1 _ 1,
, _ 1 1 _ 1
У,ц - /(У(ц), У(X(ц) - т), У(X(ц) - т), X(ц), X(X(ц) - т), X(ц)) _ 0, (2.6)
ОуУ,ц + ОхХ, ц + °/,ц _ 0 Здесь приняты следующие обозначения:
_ /у _ /х, ^ _ О _ ¿О О _ дО О =дО
У,'ц ' Х'ц /ц' Гц ф' ^ ду' °х дх' дГ
Кривая (2.1), являясь решением задачи (1.1), (1.2), в то же время является интегральной кривой системы уравнений (2.6). Для нахождения этой кривой систему (2.6) следует разрешить относительно производныхуиц, ху,ц, Хц. Успех этой операции зависит от обусловленности системы. Обусловленность же зависит от выбора параметра ц, который определяется вектором а.
В качестве меры обусловленности системы линейных уравнений примем (см. [7]) величину определителя этой системы А, деленную на произведение квадратичных норм его строк Обозначим меру обусловленности через В. В силу неравенства Адамара, для определителей абсолютная величина В удовлетворяет условию 0 < |В| < 1. Значение В = 0 соответствует вырожденной матрице системы, а |В| = 1 означает, что система уравнений является, в терминологии [7], идеально обусловленной.
Назовем наилучшим такой аргумент задачи, который доставляет системе уравнений (2.6) наилучшую обусловленность, т.е. наибольшее абсолютное значение меры обусловленности В.
Раскрыв определитель системы А по элементам его первой строки, получим
n + m +1
А = £ (-1 )k . (2.7)
k = 1
Здесь Ак - определитель, получающийся при вычеркивании в матрице системы последних п + т уравнений к-го столбца.
Тогда В = А//, где / = + т + 1 > 0. Здесь
n + m +1
«1 - ' " а2
d 1 = ^ £ а = 1,
ч+1 = Л+7Р, в = m, (2.8)
dv + n + 1 = £ gV, y + £ gV, + GV, „ v = 1, m.
Ъ=1 j=1
n
m
Заметим, что так как а - единичный вектор, то в равенствах (2.8) имеем d1 = 1, а величина d
не зависит от компонент а/ (/ = 1, п + m +1).
Согласно правилу Крамера, решение системы (2.6) можно представить следующим образом:
г +1 Аг _ ( ,4« +1 + 1 Ап + ]
= (-1 = (-1) ац у ' А' ац ^ ;
ат(-')
п + т Ап + т +1
г = 1, п, 1 = 1, т. (2.9)
(2.10)
А ' d ц А
Начальные условия примут вид
У(ц) = ф(ц), У(ц) = ¥(ц), х(ц) = v(ц), це [-цт, 0),
у( 0) = ф( Го) = Уо, х (0) = v( Г о) = хо, Г( 0) = ¿о,
где величина цт находится из решения уравнения Г0 - т = Г(-цт).
В условиях (2.10) переменная ц отсчитывается от начальной точки задачи (1.1), (1.2).
Теорема. Пусть начальная задача (1.1), (1.2) для системы дифференциально-алгебраических уравнений с запаздывающим аргументом имеет единственную гладкую интегральную кривую. Для того чтобы эту задачу преобразовать к наилучшему аргументу, необходимо и достаточно выбрать в качестве такового длину дуги А, отсчитываемую от начальной точки вдоль интегральной кривой. Тогда задача (1.1), (1.2) преобразуется к виду (2.9), (2.10) при ц = А.
Доказательство. Необходимость. Для нахождения экстремума функции В при условии (2.5) составим функцию Лагранжа
Ь =
1 X (-1)" +1 а/А/ + У 1- X
/ п+ т+1 \ 2
а/
/ = 1
/ = 1
Здесь у - множитель Лагранжа.
Необходимые условия экстремума функции Ь задаются соотношениями
дЬ
1А
= (-1 )к -2уа/ = 0, / = 1, п + т + 1, да/ а
которые удовлетворяются при а/ = (-1)/ + 1^^.
Из условия а2 = 1 определяется множитель Лагранжа
Сп + т +1 л 1/2
2
X А2
т = ±
V г =1
2 а
Таким образом, экстремум функции Лагранжа достигается при
а +1 А/
а/ = ±(-1Г
Сп + т +1 \1/2' 2
/ = 1, п + т +1.
X а2
г=1
(2.11)
Если выражения для а/ подставить в равенство (2.7), то определитель системы будет удовлетворять равенству
Сп + т +1 Л 1/2
2
X А2
А = ±
V г =1
и экстремум функции Лагранжа достигается при значениях
/ Л\/ +1 А
а/ = (-1) -т.
(2.12)
(2.13)
Сравнивая эти значения с правыми частями системы (2.9), видим, что имеют место равенства
а = Ущ ; = 1 п,
а} + п = Х],ц, ] = 1 М (2.14)
ап + м +1 = Х,ц,
которые показывают, что направление, в
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.