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

Вычислительные технологии

Том 15, № 5, 2010

Вариационные подходы к численному решению дифференциально-алгебраических уравнений

М.В. Булатов

Институт динамики систем, и теории управления СО РАН, Иркутск, Россия

e-mail: [email protected]

B.K. Горбунов, Ю.В. Мартыненко Ульяновский государственный университет, Россия e-mail: vkgorbunov@mail. ru, mar j 2005@yandex. ru

I Ii vkii Дин Конг

Институт математики Вьетнамской Академии наук и технологий, Ханой

e-mail: ndcongSmath.ас.vn

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

Ключевые слова: дифференциально-алгебраические уравнения, вариационные методы, сплайны, коллокация.

Введение

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

где х(Ь) € Яп и ^ : О0 ^ О0 С Я2п+1, 01 С Яп, Эти проблемы возникают в теории электрических цепей, механике (кинематические уравнения с голономными связями), ядерной энергетике, оптимальном управлении, а также в биологии и других областях естествознания [1]. Для таких систем, как и для ОДУ в нормальной форме

F(x'(t),x(t),t) = 0, t е [to,T]

(1)

x'(t) = f (x(t),t),

(2)

также естественно накладывать начальные условия

x(to) = x0

(3)

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

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

dF(x'(t),x(t),t)

-м--(4)

на решениях (1), В последнем случае система НОДУ называется дифференциально-ал-гебраичеекими уравнениями (ДАУ), Постановка начальных и краевых задач для таких уравнений требует учета структуры вырождения системы (1), В данной работе мы ограничимся линейными ДАУ вида

F(x'(t),x(t),t) = A(t)x'(t) + B(t)x(t) — f (t) = 0, t E [0,1], (5)

где A(t), B(t) — (n x п)-матрпцы, f (t) и x(t) — заданная и искомая n-мерные вектор-функции.

Определение [2]. ДАУ (5) имеет решение типа Коши индекса r, если существует линейная комбинация

t r

x(t,c) = Ф(1)о + i Ko(t,T )f (t )dT + £ Kj (t)f (j-1)(t), (6)

0 j=1

где $(t), K0(t,T), Kj — (n x п)-матрицы, c E Rn, гапкФ^) = R = const, V t E [0,1], обращающая (5) в тождество, и на, любом, подотрезке [а, в] С [0,1] нет решений, отличных от (6).

Дифференциально-алгебраические уравнения индекса два и выше принято называть ДАУ высокого индекса. Для таких задач ряд разностных схем, хорошо зарекомендовавших себя при численном решении жестких ОДУ, могут порождать неустойчивые процессы или быть принципиально неприменимыми в силу вырожденности матричного пучка XA(t) + B(t). Для иллюстрации вышесказанного приведем пример. Пример 1. Система ДАУ

q(t) \

, t E [0,1], а =1,

9(t) J

имеет решение типа Коши индекса два: x2 = (g' — q)/(1 — а), x1 = g — tx2. Для нее неявный метод Эйлера A(ti+1)(xi+1 — xi) + hB(ti+1)xi+1 = hf (ti+1) при —1 < а < 1 является неустойчивым, а при а = 0 неприменим из-за вырожденности матрицы A(ti+1) + hB(ti+1).

К настоящему времени разработка численных методов решения начальной задачи для линейных ДАУ идет по следующим направлениям,

1, Перенос результатов, полученных для численного решения жестких ОДУ, на ДАУ (см., например, [3-5] и приведенную там библиографию). Основной недостаток данного подхода состоит в том, что такие разностные схемы можно применять для весьма узкого

класса задач. Построены достаточно простые примеры линейных ДАУ индекса два, для которых эти методы порождают неустойчивые процессы, или принципиально неприменимы,

2, Алгебро-аналитичеекое направление (см., например, [3, 6, 7]), основанное на применении различных проекторов с дальнейшим дифференцированием для некоторых классов ДАУ, Эти алгоритмы весьма сложны в реализации и, кроме того, их можно применять только в том случае, когда у системы (4) ранг матрицы A(t) является постоянным и ее индекс не выше двух,

3, Достаточно универсальные алгоритмы: метод нормальных сплайнов (НС) [8, 9] и метод параметризации — вариационных сплайнов, которые свободны от перечисленных недостатков и применимы к более широкому классу задач. Полную информацию о таких подходах и библиографию к ним можно найти в [10-14],

Основной недостаток методов данных алгоритмов — их высокая трудоемкость, В этом плане они проигрывают традиционным конечно-разностным методам в области применимости последних. Возникает вопрос о возможности создания более экономичных методов смешанного типа1, применимых если не для произвольно вырожденных ДАУ, то для более широких классов, чем известные разностные методы. Авторы настоящей статьи, развивая идею этого направления, предлагают коллокационно-вариацион-ный метод, объединяющий идеи методов сплайн-коллокации [15] и НС [10],

1. Метод вариационных сплайнов

В данном разделе описан метод вариационных сплайнов (ВС) [9, 13] для задачи (3), (5), для которой он будет модифицирован в следующем пункте. Приближенное решение находится в виде п-мерного сплайна с подвижными узлами.

Зададим па отрезке [0, 1] сетку 0 = ¿0 <¿1 < ... < ¿м = 1 и на каждом промежутке [¿£-1, ¿о) к = 1, 2, ..., N компоненты решения будем аппроксимировать полиномами степени р:

р

Жг(*) ~ ЗД) = 7 (* - ¿—У , г = 1, . . . , П. (7)

.7=0

Далее используем векторную запись коэффициентов сплайна с£ = {.,..., с0п}. На кусочный п-мерный полином (сплайн) (7) накладываются условия непрерывности Б0 (¿о) = £>0+1^0), из которых следуют равенства

с0к+1 = £ скно, (8)

7=0

где Н0 = ¿о — ¿о-1 ■ Можно также наложить условия непрерывной дифференцируемости £0(¿0) = (¿0), из которых следуют равенства

с?+1 = £ 7-1. (9)

. =1

Из начального условия (3) следует равенство с0 = х0.

1 Вариант метода параметризации сочетает оптимизацию невязки с интегрированием регулярной управляемой системы ОДУ быстрыми одношаговыми методами.

Таким образом, сплайн (7) определяется набором параметров — векторов {с^ } и скаляров гк, которые находятся из условия минимизации интеграла невязки решаемой системы (5)

1

з (х(г)) = у ||Г (х'}х}г)ЦП ¿г, (ю)

0

где || • ||п — норма евклидова пространетва Еп, по искомым параметрам и при условии (3), Этим исходная задача (3), (5) заменена классической задачей вариационного исчисления на параметризованном классе функций (7), В случае разрешимости задачи (3), (5) минимум функционала (10) равен нулю. Соответствующую минималь будем называть "вариационным сплайном". Такие сплайны по построению отличаются от известных коллокационных сплайнов [15], аппроксимирующих решения регулярных ОДУ, Для применения метода ВС можно не делать предположения о существовании и единственности решения исходной задачи (существовании решения типа Коши индекса г), а ограничиться более слабым предположением о существовании минимума функционала (10) на выбранном классе функций. При этом функция, доставляющая положительный минимум функционалу (10), будет псевдорешением исходной задачи, В случае неединственности решения или псевдорешения в [16] указан метод регуляризации решаемой задачи, основанный на использовании пробного решения и обеспечивающий аппроксимацию решения или псевдорешения, ближайшего (в норме Ь2) к пробному решению. Мерой аппроксимации решения задачи (3), (5) будем считать значение функционала невязки (10),

Представим минимизируемый функционал (10) в виде

N ^

з(х(г)) = ^ / ||Г(х',хд)цП¿г (11)

к=1, ^

¿к-1

и будем проводить пошаговую минимизацию, требуя, чтобы уровень суммарной невязки (11) не превышал заданной величины е. Для этого достаточно обеспечить выполнение неравенств

¿к

I ||Г(х',х,г)Ц2п¿г < Нке. (12)

¿к-1

Интегральный функционал (10) и его составляющие в (11) на классе сплайнов (7) представляют собой функции конечного числа переменных {с^ ,гк}, которые при сделанных предположениях гладкости являются непрерывно дифференцируемыми. Введем, используя вид (5) функции Г (х', х, г), для каждого интервала (гк-1,гк] функцию параметров

¿к

^{скдк) = ±- I |\А(г)8'к(г) + в(г)8к(г)-№\\2пЯ- (13)

¿к-1

Назовем эту функцию частичной невязкой рассматриваемой задачи на интервале [гк-1, гк], Построение сплайна, обеспечивающего заданный уровень е суммарной невязки (11), в силу (12) сводится к последовательной минимизации функций (13) с точностью е:

<Рк (ск, г^ < е.

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

В задаче минимазации (13) число искомых параметров равно п (р + 1), однако из условия непрерывности сплайна (8) и, возможно, непрерывной дифференцируемости (9) число искомых параметров сокращается до пр и до п (р — 1) соответственно. Функция (13) с учетом (7) имеет вид

1

¿к

х

¿к-1

(ск, гк) = —х

Я-к

£ (* — *к-1 г1 (¿¿(*) + (* — *к-1)В(*)) ск + В(*)с§ — /(¿)

7=1

(14)

Дифференцируя (ск, ¿к) по ск и приравнивая частные производные к пулю, получим систему линейных алгебраических уравнений (СЛАУ)

мк Ск =

(15)

здесь Ск

к1 к1 ск , с2

т

ск ) , Мк — блочно-симметричная матрица вида

( т11 тк

12

т^ \

к к т21 т22

т:

2р

V тр1

т

К 'р2

Чр)

где тк- — (п х п)-матрицы, определенные по правилу

¿к

тк = у (* — ¿к-1)г+7-2 (¿¿(¿) + (* — ¿к-1)В(^))Т {¿¿(¿) + (* — ¿к-1)В(*))

tk-l

(16)

а вектор Ок = (^т, (т, ..., (Т) имеет компоненты ¿к

= / (* — ¿к(¿¿(¿Ж* — ¿к-1)В(^))Т (/(¿) — В(^)ск) ^ = 1,...,р. (17)

¿к-1

Реализация метода ВС требует вычисления определенных интегралов (16) и (17), Эти интегралы можно вычислять по известным квадратурным формулам, вводя на интервале интегрирована [¿к-1, ¿к] дополнительную подеетку, например равномерную tk = ¿к + 1т, I = 0,1,..., д, т = Нк /д. В приведенных численных примерах определенные интегралы вычисляются по квадратурной формуле Симпсона,

Опыт применения методов математического программирования и ВС [11-14, 16] подтвердил их эффективность при решении различных ДАУ с произвольным вырождением главной части, для которых не известно стандартных подходов, не использующих специфики решаемых задач. Однако эти универсальные методы, как и следовало ожидать, проигрывают по времени решения традиционным конечно-разностным

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

2. Метод коллокационно-вариационных сплайнов

Для ДАУ высокого индекса главная проблема построения численных алгоритмов состоит: 1 — в аппроксимации исходной задачи, 2 — в устойчивости алгоритма. Для достижения этих целей, как и в методе ВС, приближенное решение задачи (3), (5) будем искать в виде непрерывного сплайна (7), (8), Однако коэффициенты {с\,...,ск} для каждого полуинтервала \Ьк-1, Ьк) будем вычислять из условий коллокации подобно методу сплайн-коллокации [15] на некоторой подеетке Ьк-1 < Ьк1 < ... < Ьы < Ьк, где I <р:

которые при к > г — индекса исходной задачи и подходящем выборе узлов Ьк1, Ьк2, Ьы обеспечивают аппроксимацию исходной задачи. Данные условия представляют недо-определенную СЛАУ из п1 уравнений относительно пр перемениых {с^,...,ср^ : г = 1,... ,п}. Остальные переменные будем находить из условия минимума квадрата нормы сплайна в гильбертово-еоболевеком пространстве п-мерных вектор-функций ШРп[Ьк-1,Ьк]

Это условие обеспечивает устойчивость алгоритма.

Поставленная для каждого промежутка \Ьк-1, Ьк) задача минимизации (19) при условиях (18) является классической задачей квадратичного программирования (КП) относительно коэффициентов сплайна {с... : г = 1,..., п}. Конкретный вид этой задачи и свойства нового метода решения исходной дифференциальной задачи определяются следующими факторами:

р

б) выбор числа точек коллокации /;

в) расположение точек коллокации , ] = 1, 2, ..., I.

Назовем описанный метод методом коллокационно-вариационных сплайнов (КВС), В зависимости от выбранных вида сплайна и условий коллокации получим различные варианты метода КВС, Эксперименты над ДАУ индекса не выше двух и сплайнами степени не выше трех показали, что можно ограничиться равномерными коллокацион-ными сетками = Ьк + тк^ тк = Нк1, ] = 1,..., I, г < I < р, где г — индекс ДАУ (см, формулу (6)),

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

А(ьк] )вк (ьк]) + В{Ьк] )вк (ьк]) = ! (ьк]), 3 = 1,...,l,

(18)

(19)

р

п

п

Учитывая, что коэффициенты с0д находятся из условий непрерывности сплайна (8) и множитель к не влияет на аргумент минимума функции (20), переходим к более простому эквивалентному функционалу

£ о!)2114111 •

7=1

Таким образом, на каждом шаге интегрирования [¿к-1, ¿к) получаем классическую задачу КП — минимизировать квадратичную функцию (21) при ограничениях-равенствах (18), где = ¿к-1 + то, о = 1, • • •, /, тк = кк//• Методом Лагранжа она сводится к СЛАУ.

Остановимся на выборе числа точек коллокации I. Из формулы (6) видно, что решение рассматриваемой задачи содержит производные входных данных вплоть до (г — 1)-го порядка. Учитывая этот факт, число точек коллокации выбираем не менее г, Если точки коллокации выбраны равномерно, то аппроксимируем не только саму

г—1

стнгается за счет условия минимизации функции (21). Например, для параболического сплайна используется одна точка коллокации ¿к1 = ¿к, При этом получим задачу КП

II к 112 . А II к 112

НС1II +4 ||с21| ^ шт, ¿(¿к) (ск + 2с2кк) + ) (с* + сккк + сккк) = /Цк)

к, , ^2А _ ^ ^ (22)

Для кубического сплайна можно выбирать одну или две точки коллокации. В случае одной точки коллокации ¿к1 = ¿к имеем задачу

12 + 4 11 с2||2 + 36 Нс^2 ^ ш1п,

Ак (ск + 2с2кк + 3е*ку + ) (с§ + сккк + сккк + сккк) = /(¿к)•

(23)

В случае двух точек коллокации ¿к1 = ¿к-1 + кк/2 и ¿к2 = ¿к задача КП принимает вид

11|2 + 4 ||с21|2 + 36 11ск||2 ^ ш1п,

с

Ак_1/2 с\ + сккк + -скХ + Вк_1/2 ск + -с^к + -ск1г2 + = Д_1/2 , (24)

+ Дь-1/2 ^ск + ^кк + ^ I Ак (ск + 2сккк + Зсккк) + Вк (с* + сккк + скк| + сккк) = /к,

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

х' = Ах, х(0) = х, (25)

где А - скаляр, Ее А < 0, |А| > 0 [1]. Любой одпошаговый метод, примененный к этому уравнению, можно записать в виде = где г = Ак, а Я(г) — дробно-

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

Если использовать метод (22) к задаче (25), то легко заметить, что функция перехода от шага к шагу (функция устойчивости) зависит не только от г, но и от к:

_ 4(1-г) + к2(4-2г) Хг+1 ~ 4(1- г)2 + к2(г-2) Х%'

То же можно сказать и о методе (23), примененном к задаче (25), Функция устойчивости в силу ее громоздкости не приводится,

3. Численные примеры

Приведем два примера, демонстрирующих возможности методов ВС и КВС, Пример 1 [7]. Рассмотрим систему ДАУ (5), где

/ г2 2г2 г2 г ^ ( з зг г 0 ( 4г2 - зг -з

г2 г2 2г2 г , в(г) = г 3 зг 0 , I (г) = 4г2 - зг -з

г2 г2 г2 г г г 3 0 зг2 - г -з

0 0 0 0 V г г г 1 \ -зг - з

Здесь функция f (¿) соответствует точному решению х(г) = (в*, в*, в*, в*) начальной задачи с условием х(0) = (1,1,1,1). Однако на любом промежутке, содержащем точку г = 0, матрица А(г) имеет переменное вырождение, и система (5) не может быть сведена к нормальной форме, т, е, не имеет конечного индекса. Соответственно начальная задача на отрезке [0, 1] не может быть решена известными методами, кроме методов НС и ВС, В работе [7] начальная задача решалась на интервале [1, 2] (с условием х(1) = (в, в, в, в)), где система (5) имеет индекс 4, специальным коллокационно-разноетным методом. При этом на сетке с 60 узлами достигалась точность порядка 10-3, В [12] данная задача решалась методом НС, и такая же точность в пространстве Ж23 4[1, 2] достигалась при 15 узлах.

На интервале [0, 1] этот пример успешно решался методами НС [12] и ВС [11], В табл. 1 приведены результаты решения задачи методом ВС на отрезках [0, 1] и [1, 2] с уровнем невязки е = 10-8, Здесь и далее использованы следующие обозначения: ктп и ктах — минимальная и максимальная величина интервала [гк-1,гк), N — количество полученных интервалов, 7 — вычисленное значение функционала (10), А — максимальное уклонение координат полученного приближения от точного решения на построенной сетке.

Таблица 1. Пример 1

Интервал Ъ ■ 1 ¿Пип Ъ 1 'тах N 3 ас

[од] 0.1 0.2 6 7Е-9 2.6Е-3

[1,2] 0.025 0.1 16 6Е-9 7.7Е-3

Таблица 2. Пример 2. й = 0 а = 0 д = 1

к 0.1 0.05 0.025

а22 1.3Е-1 7.1Е-2 3.7Е-2

А2з 1.3Е-1 7.1Е-2 3.7Е-2

А24 8.0Е-4 2.1Е-4 5.6Е-5

Таблица 3. Пример 2. й = 0 а = 0 д = 0

к 0.1 0.05 0.025

а25 * * *

а26 * * *

а27 1.2Е-3 3.4Е-4 5.7Е-5

Метод КВС был реализован для параболических и кубических сплайнов. Для них основная для метода задача КП имеет варианты (22), (23) и (24), Эти программы были протестированы на ряде примеров, один из которых представлен далее. Пример 2. Рассмотрим систему ДАУ

( 1 Ь N х + ( 0 Я \ х = ( ехр(Ь) + (д — ¿)ехр(—Ь) N ^ 0 + ^ 1 Ь + а/х ^ехр(*) + (* + а — ^)ехр(—Ь) )

на отрезке Ь € [0,1] с начальным уеловием х1(0) = х2(0) = 1, соответствующим решению х = (ехр(Ь), ехр(—Ь))т, Индекс системы зависит от значений параметров а, д. Если d = 0 т0 матрица

обратима и система имеет индекс, равный нулю, т, е, сводится к ОДУ без дифференцирования, Нетрудно проверить, что при d = 0 и а = 0 индекс равен единице и при d = 0, а = 0 q = 1 _ двум. Результаты расчетов для последнего набора параметров методом КВС представлены в табл. 2, Здесь и далее Д22, Д2з и Д24 — погрешности методов (22), (23) и (24).

Эта же задача решалась методом ВС. При выборе сплайна степени 3 и уровня невязки е = 10-11 сформировалась сетка N = 16 hmin = 0.05 и hmax = 0.01. Погрешность аппроксимации первой координаты на удвоенной сетке равна 0.38, второй координаты — 0.23. Отметим, что здесь точность заметно хуже, чем в методе КВС. d = 0 а = 0 q = 0

ричный пучок AA(t) + B(t) сингулярен. В этом случае применение ряда неявных многошаговых методов, наиболее эффективных для решения ДАУ, принципиально невозможно. Результаты расчетов такой задачи методом КВС представлены в табл. 3. Знак * означает, что метод КВС с данным вариантом примера дает очень большие погрешности.

Заключение

Представленные вариационные методы ВС и КВС, как и упомянутые во введении методы НС и метод параметризации, имеют очевидные преимущества над ранее разработанными разностными методами (см. работы [12, 15] и библиографию к ним). Их можно применять для численного решения более широкого класса задач, в том числе для ДАУ, содержащих особые точки. Приведенные примеры демонстрируют универсальность (относительно классов решаемых задач) метода ВС. Однако он достаточно затратный и уступает методу КВС при усложнении структуры вырождения главной части (в рамках систем конечного индекса, где можно применять метод КВС). Таким образом, оба метода — ВС и КВС — имеют не перекрываемые преимущества и перспективны для дальнейшего развития.

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

полагать, что возникающие при этом СЛАУ будут иметь числа обусловленности, не зависящие от степени сплайна p, а определяемые только шагом сетки и структурой самой задачи,

Нгуен Дин Конг благодарит Вьетнамский Национальный Фонд Науки и Развития Технологий (NAFOSTED) за поддержку данных исследований.

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

[1] Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи. М.: Мир, 1999. 685 с.

[2] Булатов М.В., Чистяков В.Ф. Применение коллокационных методов для решения сингулярных линейных систем ОДУ // Модели и методы исследования операций. Новосибирск: Наука, 1988. С. 164-170.

[31 Булатов М.В. Методы решения дифференциально-алгебраических и вырожденных интегральных систем: Дисс. ... д. ф.-м. н. Иркутск: ИДСТУ СО РАН, 2002. 244 с.

[4] Brenan К.Е., Campbell S.L., Petzold L.R. Numerical Solution of Initial-Value Problems in Differential-Algebraic equations. North Holland. New York, 1989. 256 p.

[5] Булатов M.В., Чистяков В.Ф. Об одном численном методе решения дифференциально-алгебраических уравнений // Журн. вычисл. математики и мат. физики. 2002. Т. 42, № 4. С. 459-470.

[6] Чистяков В.Ф. Системы интегро-дифференциальных уравнений с тождественно вы-

...

[7] Бояриицев Ю.Е., Данилов В.А., Логинов A.A., Чистяков В.Ф. Численные методы решения сингулярных систем. Новосибирск: Наука, 1989. 223 с.

[8] Горбунов В.К. Метод параметризации задач оптимального управления // ЖВМ и МФ. 1979. Т. 19. № 2. С. 292-303.

[9] Горбунов В.К., Мартыненко Ю.В. Метод вариационных сплайнов для неявных дифференциальных уравнений // Вестник Самарского госуд. технич. ун-та: Серия "Математическая". 2007. № 2(6) С. 16-28.

[10] Горбунов В.К. Метод нормальной сплайн-коллокации // Там же. 1989. Т. 29, № 2. С. 212-224.

[11] Горбунов В.К., Лутошкин И.В. Развитие и опыт применения метода параметризации в вырожденных задачах динамической оптимизации // Изв. РАН. Теория и системы управления. 2004. № 5. С. 67-84.

[12] Горбунов В.К., Петрищев В.В. Развитие метода нормальной сплайн-коллокации для линейных дифференциальных уравнений. // Журн. вычисл. математики и мат. физики. 2003. Т. 43, № 8. С. 1161-1170.

[13] Gorbunov V.K., Lutoshkin I.V. The parametrization method in optimal control problems and differential-algebraic equations //J. Comput. Appl. Mathem. 2006. Vol. 185, is. 2. P. 377-390.

[14] Gorbunov V.K., Sviridov V.Yu. The method of normal splines for linear DAEs on the number semi-axis // Appl. Numer. Math. 2009. Vol. 59, is. 3-4. P. 656-670.

[15] Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.

[16] Gorbunov V.K., Lutoshkin I.V., Martynenko Y.V. A parametrization method for the numerical solution of singular differential equations // Appl. Numer. Math. 2009. Vol. 59, is. 3-4. P. 639-655.

Поступила в редакцию 12 октября 2009 г., с доработки — 20 апреля 2010 г.