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

ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА

УДК 519.62

doi: 10.18101/2304-5728-2017-3-3-9

О КОЛЛОКАЦИОННО-ВАРИАЦИОННОЙ РАЗНОСТНОЙ СХЕМЕ ДЛЯ ДИФФЕРЕНЦИАЛЬНО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ1

© Соловарова Любовь Степановна

кандидат физико-математических наук,

младший научный сотрудник,

Институт динамики систем и теории управления

имени В. М. Матросова СО РАН

Россия, 664033, Иркутск, ул. Лермонтова, 134

E-mail: soleilu@mail.ru

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

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

Введение

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

1 Работа выполнена при финансовой поддержке РФФИ, грант 16-31-00219.

3

могут быть неустойчивыми. Другие направления теории численного решения ДАУ (так называемые «расширенные системы» [7] и подходы, основанные на различных полу обратных матрицах и проекторах [11]) трудны в реализации и не всегда применимы.

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

Определение 1. Жесткая система ОДУ — система, для которых существуют участки области определения с различными скоростями изменения компонент решения.

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

x\t) = Xx(t),x(Qi) = x0,t&[<d,\], (1)

где Л — скаляр и Re Л « 0 .

Любой одношаговый метод для задачи (1) можно записать в виде хм = R{z)xi, где z = Лк, R(z) — полином или дробно-рациональная функция, которую принято называть функцией устойчивости.

В дальнейшем нам потребуется следующее определение. Определение 2. [6] Метод называется А -устойчивым, если |i?(z)| < 1 для

любого z, лежащего в левой полуплоскости. Если lim_ /¿(z) = q < 1, то метод называется сильно А -устойчивым.

В настоящей работе для численного решения линейных ДАУ, в том числе и жестких, предложена коллокационно-вариационная разностная схема. Данная схема является одним из вариантов коллокационно-вариационного подхода [4], частные случаи которого изложены в [1], [2], где в качестве координатных функций выбраны непрерывные сплайны степени не выше трех.

1. Постановка задачи

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

A{t)x\t) + B(t)x(t) = f{t), t e [0,1], (2)

x(0) = x0, (3)

где A(t) и В(/) — (n x n) -матрицы. /(0 й x(t) — заданная и искомая «-мерные вектор-функции, причем dct А(/) = 0. Как уже отмечалось, такие задачи принято называть дифференциально-алгебраическими уравнениями. Предполагается, что для ДАУ (2) начальное условие (3) согласовано с правой частью, т. е. рассматриваемая задача имеет достаточно гладкое решение.

Характеристикой сложности ДАУ является понятие индекса, которое имеет несколько определений. В случае вещественно-аналитических коэффициентов матриц A(t), B(t) эти определения эквивалентны [8].

Определение 3. [8] Система (2) имеет семейство решений типа Коши

Жс) = Ф(Ос + }го(^т)/(т)£/т + Х^ДО/ОЧ)(0, (4)

О 3= о

где Ф(0 , K0(t,z) ,Kx{t) ,...,Kr(t) — (п х п) -матрицы, rank<S(t) =

= г = const We[0,l],cei?" и на любом интервале [a,/?]ç[0,l] нет решений, отличных от (4). В этом случае г принято называть индексом исходной системы (2).

Если г > 2, то такие задачи принято называть ДАУ высокого индекса. Разработка методов численного решения таких задач и жестких ДАУ индекса один наталкивается на большие трудности. Для иллюстрации вышесказанного приведем примеры. Пример 1.[ 9]

ОМ'ЛГло'

,1 at J[v(i)J U(0/

где а — числовой параметр. Данный пример имеет единственное решение u(t) = g(t)~ atv , v(t) = f(t)~ g'(t). Индекс данного ДАУ равен 2.

Зададим на отрезке интегрирования [0,1] сетку /( = ih. i = 0.1..... N . N = 1/ h и обозначим Д = A(tt), fi = B(r), fr= f(tt), x, « x(tt) .

Применим для данного примера неявный метод Эйлера, который имеет вид

Ам(хм -x,) + hBt+XX1+1 =|/Л ,-где I// , =h(f(ti+l),g(ti+l))T, хм =(ui+l,vM)T *(u(tM),v(ti+l)), u0=f{V)-g'{t)\t=0,u0=g{Q). Опуская выкладки, получим

(\ aiVM'(0Vf0 l + aY«(i) 0 0 Jl v'(0

а

1

u,+i=g,+i-af,+iv,+i>v,+i=- , ,

1 + a 1 + a

При a < -0.5 будем иметь неустойчивый процесс. Пример 2. [10]

(\ -at\(u'{t)\_( к a(\-Xt)Yu(t)^ 0 0 llv'(0j"l-l 1 + cct llv(0

Ге[0,1],

где а и Л — скалярные параметры. Точное решение данного примера м(0 = (1 + а0ехр(/1/) = ехр(Л^ . При Л « 0 задача жесткая. Опуская выкладки, получим, что неявный метод Эйлера для данного примера дает рекуррентное соотношение

( и ^

= R(z,a>)

f 0 l + at.AfuA

о

1

v.

/V 1 у

где R(z,co) = (1 -со) / (1 - z -со), со = ah, z = ЛЬ . Видно, что для любого

г <0 можно подобрать со . что > 1, т. е. метод будет неустойчи-

вым. В статье [10] показано, что для численного решения данного примера ряд неявных методов Рунге — Кутты и многошаговых алгоритмов также неэффективен: для устойчивости требуется достаточно малый шаг интегрирования.

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

2. Коллокационно-вариационная разностная схема

Данный параграф посвящен построению коллокационно-вариационнной разностной схемы с одной точкой коллокации.

Пусть Ь2(хм,х^х1ч,0 ,1 — интерполяционный вектор-

полином, проходящий через точки , ч), , ), , ) • Подставив его в (2) и полагая / = I/ ,, получим

Ам (Зх,+1 - 4х, + хн) + 2кВмхм = 2к/м. (5)

Будем смотреть на разностную схему (5) как на одношаговую, т. е. будем иметь недоопределенную СЛАУ размером п х 2п , которую дополним условием минимума выпуклой целевой функции Ф(х1+1,х1).

Определим ||Х2 (хм, х{ хя, , ^ е |/, ,, ] как

2 ',+1

\\12 0|| = Е I 12]]Т ' х<' ' 04') 5 Х, 5 Х,-1 > т. (6)

Данная функция аппроксимирует квадрат нормы решения на отрезке

\fi-M ■

Применяя для вычисления определенного интеграла (6) формулу левых прямоугольников, будем иметь

2 <«

(') = X, (Х1+1 > Х1' Х1-1' 04 (Х1+15 Хг ■> Х1-1' 0^ '

3=0 I

(7)

■■2к

—(-х1+1 +4х,-3х,,1) 2 п

2

—(х,+1 -2ж,.+*,._!)

Здесь норма вектора — евклидова.

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

ф(х,+1 ,*,■)=—|-х,+1+4х, - 3х,-11Г+1х,+1 -2х,+х,-11Г •

(8)

Выражение (8) представляет собой некий аналог квадрата нормы интерполяционного векторного многочлена Ь2 (х1+1, xi, хя, 0 , (е |/, ,, ti+l ].

Стандартно [3] перейдем от задачи на минимум функции (8) с ограничением (5) к задаче на безусловный минимум функции Лагранжа

г(хм,х{,А) = —|-х,+1 + 4х, -Зх^Ц2 +|х,+1 -2х, +Х,_1|2 +

+АТ[(ЗАМ +2кВм)хм -4Д+1х, + А^-2к/м], где А = (Л1,Л2,...,Лп)г — «-мерный вектор множителей Лагранжа. Аргументы х1+1,х15Л, при которых функция (8) будет достигать минимума, удовлетворяют следующей СЛАУ

((2 + 0.5 к2)Е -(4 + 2 к2)Е {ЪА1+1+2кВм)т\х1+^ -(4 + 2 к2)Е (8 + 8/2 2)Е -4 Атм х,

ЗД.+1+2А5<+1 -4Д+1 0 1А

Г-(2- И.5/2 2)Ел ( 0 ^

(4- ±6к2)Е + 0

V ~Д+1 )

Таким образом, полагая хм известным, х/ и х/ 1 находим из условия минимума функции (8) при ограничениях (5), т. е. продвигаемся сразу на 2к . Далее процесс повторяется до исчерпания отрезка.

Опуская выкладки, приведем функции устойчивости, которые дает предлагаемая схема для задачи (1)

х, =^(г,/2) =

0.25/г (4 - Зг)(1 - г) + (1 - г)(1 - 2г)

4/г2 (г -1)2 - (2г -1)2

-4/72(г-1)-(2г-1)

4/г2 (г -1)2 - (2г -1)2

Заметим, что данные функции устойчивости отличаются от функций устойчивости всех известных методов [5], так как зависят не только от г, но и от к. Непосредственные выкладки показывают, что 1пп _ /^ (г.И) = 1 / 2 . 1пп _ = 0 . т. е. данный алгоритм в совокуп-

ности является сильно А -устойчивым.

3. Численные расчеты

В данном параграфе приведем результаты численных расчетов метода (8), (5) для модельных примеров, приведенных выше.

Для примера 1 в качестве точного решения возьмем и = ехр(/) -а/(ехр(0.5/) -ехр(0), V = ехр(0.5/) -ехр(/), тогда правая часть примет вид /(0 = ехр(0.50 , g(t) = exp(t). Расчеты для этого случая при а = -0.6 приведены в таблице 1.

Таблица 1

h 0.1 0.05 0.025 0.0125

erru 0.1 0.06 0.036 0.019

errv 0.175 0.1 0.06 0.031

Здесь и ниже егги = тах, , ,, \и1 - , еп; = тах, , ,, - .

Для примера 2 значения параметров выберем следующим образом: а = 100. /1 = 10. Результаты численных расчетов приведены в таблице 2.

Таблица 2

h 0.1 0.05 0.025 0.0125 0.00625

erru 0.183 0.14 0.09 0.049 0.025

errv 0.045 0.035 0.022 0.012 0.006

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

Заключение

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

Литература

1. Вариационные подходы к численному решению дифференциально-алгебраических уравнений / М. В. Булатов, В. К. Горбунов, Ю. В. Марты-ненко, Нгуен Дин Конг // Вычислительные технологии. 2010. Т. 15, № 5. С. 3-14.

2. Булатов М. В., Рахвалов Н. П., Соловарова Л. С. Численное решение дифференциально-алгебраических уравнений методом коллокационно-вариационных сплайнов // Журнал вычислительной математики и математической физики. 2013. Т. 53, № 3. С. 46-58.

3. Васильев Ф. П. Численные методы решения экстремальных задач. М.: Наука, 1988. 552 с.

4. Горбунов В. К. Метод нормальной сплайн-коллокации // Журнал вычислительной математики и математической физики. 1989. Т. 29, № 1. С.212-224.

5. Деккер К., Вервер Я. Устойчивость методов Рунге — Кутты для жестких нелинейных дифференциальных уравнений. М.: Мир, 1988.

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

7. Чистяков В. Ф. О расширении линейных систем, не разрешенных

относительно производных. Иркутск, 1986. 25 с. (Препринт ИрВЦ СО АН СССР;5).

8. Чистяков В. Ф. Алгебро-дифференциальные операторы с конечномерным ядром. Новосибирск: Наука, 1996. 280 с.

9. Brenan К. F., Campbell S. L., Petzold L. R. Numerical solution of Initial-Value Problems in Differental-Algebraic Equations. Philadelphia: Appl. Math., 1996.270 p.

10. Kunkel P., Mehrmann V. Stability properties of differential-algebraic equations and spin-stabilized diskretizations // Electr. Trans. Numer. Analys. 2007. Vol. 26. P. 385-420.

11. Lamour R., März R., Tischendorf С. Differential-Algebraic Equations: A Projector Based Analysis. Springer, 2013.

ABOUT COLLOCATION-VARIATIONAL DIFFERENCE SCHEME FOR DIFFERENTIAL-ALGEBRAIC EQUATIONS

Lubov S. Solovarova

Cand. Sei. (Phys. and Math.), Junior Researcher,

Matrosov Institute for System Dynamics and Control Theory SB RAS

134 Lermontova St., Irkutsk 664033, Russia

The article deals with the systems of ordinary differential equations with an identically degenerate matrix preceding a derivative (differential-algebraic equations). The initial condition is supposed to be consistent with the right-hand side, so the problem under consideration has a solution. We have given necessary definitions from the theory of solving differential-algebraic and stiff ordinary differential equations and shown the difficulties of their numerical solution. A special case of the collocation-variational method for solving differential-algebraic equations containing stiff components and having index not exceeding two is considered. We have obtained the stability functions for the Dahlquist model equation, and made numerical calculations of the known test examples.

Keywords: differential-algebraic equations; difference schemes; index; stiff problems; A-stability; stability functions.