CC BY
32
Вычислительные технологии
Том 10, № 2, 2005
О СИСТЕМАХ НЕ ТИПА КОШИ - КОВАЛЕВСКОЙ ИНДЕКСА (1,k)
С. В. ГАйДОМАК, В. Ф. Чистяков Институт динамики систем и теории управления СО РАН,
Иркутск, Россия e-mail: [email protected], [email protected]
An approach to the investigation of the non Cauchy — Kowalewska type systems is suggested. Systems of the first order partial differential equations with variable coefficients are considered in which a singular matrix stands in front of the evolutionary term. A classification of such systems is given, based on the representation of the matrices, which specify the system. A definition of systems having the index (1, k) is introduced. The
method of straight lines for the numerical solution of some classes of the systems is
validated.
Введение
Рассматривается начально-краевая задача вида
Л1;1м := ADtu + BDxu + Cu = f; (1)
u(x,t0) = ф(х), u(x0,t) = ^(t), (2)
где A = A(x,t), B = B(x,t), C = C(x,t) — (n x п)-матрицы с элементами, зависящими от переменных x и t; (x,t) £ U = [x0,X] x [t0,T] С R2; u = u(x,t) — искомая, f = f(x,t), ф^), ^(t) — заданные вектор-функции; Dt = d/dt, Dx = d/dx. Предполагается, что входные данные задачи (1),(2) достаточно гладкие и
det A(x,t) = 0 V(x,t) £ U. (3)
Допускается одновременное вырождение матриц A и B и, более того, вырождение пучка матриц Л A + B для любых (x, t) £ U и Л £ R.
Система (1), удовлетворяющая условию (3), называется системой не типа Коши — Ковалевской [1]. Системы не типа Коши — Ковалевской возникают в различных областях приложений, в частности, они встречаются в гидродинамике, газовой динамике, теории малых колебаний жидкости, теории тепломассообмена [1-3].
В большинстве работ, посвященных системам (1), удовлетворяющим условию (3), рассматриваются уравнения с постоянными матрицами коэффициентов. Изучение систем с переменными матрицами коэффициентов связано с особыми трудностями. Авторы настоящей работы предлагают подход к исследованию систем не типа Коши — Ковалевской, основанный на использовании методов из теории алгебро-дифференциальных систем.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.
1. Основные определения и вспомогательные утверждения
Рассмотрим систему уравнений
(ДДх + В)и = С, (4)
где А = А(х, ¿), В = В(х, ¿) — (п х п)-матрицы; £ = С(х,£) — заданная п-мерная вектор-функция. Предполагается, что А, В, £ обладают требуемой при рассуждениях гладкостью.
к
Определение 1. Оператор Лк := ^ ¿¿(х,£)ЦХ такой, что
г=0
Лк о (АД^ + В)и = Д^и + Лк[В(х, ¿)]и Уи € Ск(и)
называется левым регуляризирующим оператором (ЛРО) порядка к, где Д(х,£) — некоторые непрерывные, по крайней мере, матрицы в области и. Наименьшее число к назовем индексом системы (4).
Если А, В не зависят от х, то для существования ЛРО необходимо и достаточно, чтобы старший коэффициент многочлена det(AA + В) = а^(£)А^ + ... не имел нулей на отрезке [¿о,Т]. Если матрицы А, В зависят и от х, то необходимое и достаточное условие существования ЛРО можно найти в монографии [4].
Системы (4), для которых определен ЛРО, разрешимы при любой достаточно гладкой правой части £ и имеют конечномерное многообразие решений. Точнее говоря, справедливо следующее утверждение.
Теорема 1. Если существует ЛРО для системы (4), то ее общее решение можно записать в виде
х к—1
и = Ф(х, ¿)с(£) + £(х, ¿), ^(х,£) = I К(х,£,з)((х, з)^з + '¡Те, (х.О^Х С. (5)
.7=0
хо
где Ф(x,t), K(x,t,s), Cj(x,t), j = 0,1,...,k — 1, — некоторые, по крайней мере, гладкие в области U матрицы; c(t) — произвольная вектор-функция. Если многообразие решений имеет нулевую размерность ^(x,t) = 0), то решение системы (4) единственно и
к-1
u = Y1 Cj (x,t)Dx c. (6)
j=0
Теорема 1 доказана в работах [4, 5], когда матрицы A, B и вектор-функция Z зависят только от x, а доказательство теоремы 1 получается подстановкой входных данных, зависящих от параметра t, в соответствующие формулы.
Теорему 1 используем для выделения класса систем (1). Будем предполагать, что в системе (1) rank A = d = const. Тогда найдутся такие невырожденные в области U матрицы
P = P(x,t) и Q = Q(x,t), что при умножении системы (1) слева на матрицу P и замене
переменной u = Qz система примет вид
Ed 0 \ f B11 B12 N n , f C11 C12 \ _ ^
0 0 ) Dtz + ( B21 B22 j DxZ + I C21 C22 )z g (7)
где д = Р/; В ,, С ,, г, ^ = 1, 2, — блоки соответственно матриц РВ^ и РАР4^+РВРЖ^+ РОД
Второе уравнение из системы (7) можно записать в виде
(В22Рж + С22)г2 = /2, (8)
где /2 = д2 - В21^ж^1 - С21 ^1, (^7 г^)т = г, (д^ д^)т = д.
Предположим, что для системы (8) определен ЛРО, тогда из нее можно выразить вектор-функцию г2 через г1, используя формулу общего решения (5). Результат можно подставить в первое уравнение системы (7). Получим систему интегродифференциальных уравнений, разрешенную относительно эволюционного члена:
к х Р^ + В,(х,£)РЖг1 + С(х, £, в)г1(в, £)^з = р(х, £), (9)
7=0 ^
жо
где матрица Ф(х, £) и вектор-функция £(х, £) получены по формуле (5) применительно к системе (8), р(х,£) = д1 — (В12РЖ + С12) [Ф(х,£)с(£) + £(х,£)].
Замечание 1. Имеются существенные различия при постановке задачи Коши для системы (1), когда Ф(х,£) = 0 или Ф(х,£) = 0. В первом случае уравнение (9) является дифференциальным, и поэтому для исходной системы (1) имеет смысл (по крайней мере формально) постановка задачи Коши. Во втором случае для выделения единственного решения необходимо задание начальных условий не только на оси Ох, но и на оси О£.
По аналогии с системой (4) введем следующее определение.
Определение 2. Будем говорить, что система (1) имеет индекс (1,к) : 1 — по переменной £ и к — по переменной х, если существует ЛРО порядка к для оператора
В22Рж + С22.
Системы с индексом (1,к), например, естественным образом возникают при переходе
т
от уравнения высокого порядка Р^ + ^ а, РЖV = 0 к системе первого порядка вида (1) с
7=0
искомой вектор-функцией и = (V Р^ ... Ржт—мт и коэффициентами
А = ^{1,0,..., 0}, В =( а ^ , N =( Ет—2 0) , С = ^аё{0, —Ет—1}, где Ь =(10 ... 0)т, а = (а1 а2 ... ат). Здесь ЛРО для соответствующей системы (8) имеет
т— 1
вид £ N7 РЖ+1.
7=0
Уравнение (9) без интегральной части является линейным дифференциальным уравнением порядка к, которое может быть гиперболическим, параболическим или эллиптическим. В целях корректной постановки задачи для уравнения (1) важно знать тип, в данном случае — уравнения (9). Тип системы (9) отнесем и на исходную систему (1). В дальнейшем будем рассматривать только гиперболические уравнения, поэтому необходимы признаки, указывающие на гиперболичность систем.
2. Признаки гиперболичности систем индекса (1, к)
Ниже нам потребуются следующие определения и утверждения.
Определение 3. Пучок матриц AA + B удовлетворяет критерию “ранг-степень” в
параметр из R1.
Определение 4. Если: 1) rank A = d = const, rank S = d + l = const V(x, t) G U, S =
в нуль ни в одной точке области U, то пучок матриц AA + ^B + C удовлетворяет двойному критерию “ранг-степень” в области U.
Определение 5. Матрица A- называется полуобратной к матрице A, если она удовлетворяет уравнению AA-A = A.
Замечание 2. В случае l = n — d, En — SS- = 0, где S = (A B), S- — любая полуобратная матрица к матрице S : SS-S = S, двойной критерий совпадает с критерием “ранг-степень”.
Лемма 1. Пучок матриц AA + B удовлетворяет критерию “ранг-степень”, тогда и только тогда, когда найдутся такие невырожденные в области U матрицы P1 = P1(x,t) и Q1 = Q1(x, t), что
где Е^, Еп—^ — единичные матрицы размерностии d и п — ^ соответственно; 7 = 7(х, £) — некоторая квадратная матрица размерности d.
Для случая, когда матрицы А, В зависят от одной переменной, лемма 1 доказана в [4, 6]. Для случая двух переменных доказательство аналогично.
Лемма 2. Если пучок АА+^В+С удовлетворяет двойному критерию “ранг-степень”, то найдутся такие невырожденные в области и матрицы Р2 = Р2(х,£) и ^2 = ^2(х,£), что
/Е^ 0 0 \ / J 0 B^ \ / C11 C12 0 \
P2(AA + ^B + C)Q2 = A I 0 0 01+ 0 Ei 0 I + I C21 C22 0 I , (11)
\ 0 0 0 / \ 0 0 0 / \ 0 0 ECTy
где Е^, Ег, Ест — единичные матрицы размерности d, I и а = п — d — I соответственно; 7, Вй, С,, г = 1, 2,^’ = 1, 2, 3, — некоторые матрицы переменных х,£ подходящей размерности.
Доказательство. Лемма 2 анонсирована в [7] без доказательства для случая, когда матрицы пучка АА + ^В + С зависят только от одной переменной. Поэтому приведем здесь доказательство. В силу постоянства ранга матрицы Б найдется невырожденная для всех (х, £) € и, обладающая свойством
В матрице LS число нулевых строк равно а = n — l — d. Далее, в силу определения 4 блок C2 имеет полный ранг для любых (x,t) G U, иначе говоря, в точке, где rank C2 < а, определитель det(AA+^B+C) = 0 VA, Следовательно, найдется невырожденная (nxn)-
матрица R со свойством C2R = (0 Ест). В силу постоянства rankA2 найдется матрица L1, такая, что
области U, если deg{det(AA + B)} = rank A = d = const V(x,t) G U, где A — числовой
(A B); 2) det(AA + ^B + C) = a0Ad^1 + ..., где коэффициент a0 = a0(x,t) не обращается
(10)
где число нулевых строк равно n — d. Тогда умножением слева и справа на матрицы diag{Li, ECT}L и R пучок матриц приведется к виду
/ Ац A12 Аіз\ / B11 B12 ВіД / C11 C12 С1з\
ЛІ 0 0 0 1+ ^ I B21 B22 B23 I + I C21 C22 C23 I .
\ 0 0 0 / \ 0 0 0 ) \ 0 0 E-J
Далее, определитель
j-M/ti D i m A 4- f ЛА11 + ^B11 + C11 ЛА12 + ^B12 + C12^\
det(AA + ^ + C) = q de4 MS21 + C21 MS22 + C22 )’
где q = det(diag{L1, ECT}LR). Степень блочного определителя по Л равна d по условию леммы. С другой стороны, если rank(A11 A12) < d в некоторой точке области U, то и степень по Л меньше d в этой точке. Полученное противоречие означает, что rank(A11 A12) = d V(x,t) Є U. Тогда найдется матрица і? со свойствами
detR=0, (A11 t)R=(Ed S).
После умножения блочной матрицы, стоящей под знаком определителя на R справа, мы придем к пучку ( ) ( )
Л /Ed 0А + (^B11 + C11 ^B12 + <12^
у 0 0J \^B21 + C21 ^B22 + С22/
Для того чтобы определитель последней матрицы имел степень d по Л, необходимо и
достаточно, чтобы определитель det(^B22 + C22) ф 0. При этом старший коэффициент определителя с необходимостью имеет вид det(^B22+C22)Лгі+• • • Он имеет степень, равную l по ^ V(x, t) Є U тогда и только тогда, когда det(B22) = 0 V(x, t) Є U. Итак, умножая пучок слева на матрицу diag{L1,E-}L и справа на матрицу Rdiag{R, ECT}diag{Ed,B—1, ECT}, мы приведем его к виду
/Ed 0 0\ i -¿^11 B12 B13\ / (S'il <^12 СЦ
Л I 0 0 0 I + ^ I Ü21 El B23 І + I <21 <22 <23 I •
0 0 0 0 0 0 0 0 E-
Умножая последнюю строку пучка на соответствующие множители и вычитая из первых двух, мы обратим в нуль блоки C13, C23. Аналогичными действиями преобразуем матрицу при ^ к искомому виду. Лемма доказана. □
Теорема 2. Система вида (1) имеет индекс (1, 0) тогда и только тогда, когда пучок ЛА + B удовлетворяет критерию “ранг-степень” в области U. Если при этом все корни многочлена det^A + B) вещественные и простые, то система гиперболическая.
Система вида (1) имеет индекс (1,1), если пучок матриц ЛА + ^B + C в области U удовлетворяет двойному критерию “ранг-степень”. Если при этом все корни многочлена det^A + D), D = B + (En — SS-)C, где S = (A B), S- — любая полуобратная матрица к матрице S, вещественные и простые, то система гиперболическая.
Доказательство. При приведении системы (1) к виду (7) положим P = P1 и Q = Q1, где Pi и Qi — матрицы из леммы 1. Индекс системы (8) равен нулю по определению тогда и только тогда, когда det B22 = 0 V(x,t) Є U .В условиях теоремы B22 = En-d. Корни многочлена deg{det^A + B)} вещественные и простые. Они совпадают с собственными
числами матрицы 3, взятыми с обратным знаком. Следовательно, система является гиперболической .
Для доказательства второй части теоремы положим Р = Р2 и Q = где Р2 и
^2 — матрицы из леммы 2. Индекс системы (8) равен единице тогда и только тогда, когда пучок ЛВ22 + С22 удовлетворяет критерию “ранг-степень” [6]. В условиях теоремы
В22 = ^ В23^ , С22 = и deg{det(ЛB22 + С22)} = I = гапкВ^. Здесь учте-
но, что последние а = п — (^ + /) строк матриц Р2АБ^ и Р^Д^ являются нулевыми. Следовательно, система (1) имеет индекс (1,1).
Покажем, что система (1), если корни характеристического многочлена det(ЛA + Б) вещественные и простые, гиперболическая. Умножим пучок ЛА + Б слева и справа на матрицы Р2 и Q2 соответственно и рассмотрим матричное выражение
ЛP2AQ2 + P2BQ2 + Р2(Е — ББ )CQ2,
от которого перейдем к следующему:
ЛАр2А^2 + Р2 BQ2 + (Е — Б,)P2CQ2• (12)
/^11 ^12 ^ П \
Здесь ,9 = Р2Б<32<32 1Б-Р2 1 = ^21 ^22 ^23 , <32 = ( ^ ^ ) и ГапкБ = d + /, где ^ -
V П П П) VП ^
некоторые блоки подходящей размерности. Тогда в силу [4] гапк(Е — Б) = п — ^ + /), поэтому
/П П ¿>13 \
Е — Б = І П П >23 I •
\П П Ега-(гі+і))
Следовательно, (12) с учетом условия теоремы можно записать
(Ел П П\ /3 П В1Д /П П >13 \
Л П ПП + П Е П + ПП >23 • (13)
\П П п) \п п П ) \П П Е„-(гі+0у
Из (13) следует, что все собственные числа матрицы 3 вещественные и простые. В силу того, что пучок матриц ЛА + + С удовлетворяет двойному критерию “ранг-степень”
и по лемме 2 он приводится к виду (13), характеристики уравнения (1) вещественные и простые, следовательно, система (1) гиперболическая. Теорема доказана. □
Сформулируем еще одно утверждение.
Лемма 3. Если в начально-краевой задаче
X
Ат + ЖДХ№ + Ьт + J К(х, з, ¿)т(з, = Л, т(х,£0) = ^(х), т(х0,£) = ф(£), (14)
хо
где Ж = Ж(х,і), Ь = Ь(х,і), К = К(х,з,і) — (d х d)-маmрицы, все собственные числа матрицы Ж положительные и простые в области и, начальные и краевые данные согласованы в точке (х0,і0) до второго порядка включительно, в частности ф(£0) = ^(х0), то (14) имеет решение т = т(х,і) в области и при достаточной гладкости входных данных.
Доказательство практически полностью повторяет рассуждения из [8], где обоснование проводится по схеме доказательства существования решения у симметрической t-гиперболической (по Фридрихсу) системы. Оно опирается на построение в области U разностной сетки с шагами h, т по пространственной и временной переменным таким образом, чтобы в этой области укладывалось целое число прямоугольников со сторонами h и т .В уравнении (14) вводится новая переменная v(x,t) = w(x,t) — 0(t) и после аппроксимации: Dtv на Atv = (v(x, t + т) — v(x, ¿))/т, на Axv = (v(x, t) — v(x — h, t))/h,
ж i—i
/ K (x, s, t)v(s, t)ds на hJ^K (x, jh, t)v(jh, t), где Ax = A/Ax, At = A/At, получаем си-
жо j=0
стему разностных уравнений
i—1
где g(x, t) = h(x,t) — 0'(t) — L(x,t)0(t) — K(x, jh, t)0(t)ds, которая позволяет получить
j=0
приближенные значения в узлах сетки v(jh,iT) = vij- (сеточную функцию). Проводится интерполяция сеточной функции функцией
v(x,t) = vij (i + 1 — x)(j + 1 — i) + vi/j (x — i)(j + 1 — i) + vij/ (i + 1 — x)(i — j)vi/j/ (x — i)(i — j),
где i' = i + 1, j' = j + 1, x = x/h, i = t/т. При условии h/т = const показывается, что семейство функций {v(x,t)} компактно в смысле равномерной сходимости и предел этого семейства удовлетворяет начально-краевой задаче (14).
Перейдем к теореме существования решения системы (1).
Теорема 3. Пусть для системы (1) выполнены условия: 1) пучок матриц AA+^B + C удовлетворяет двойному критерию “ранг-степень”; 2) корни многочлена det(AA + D), D = B + (En — SS—)C отрицательные и простые, S = (A B), S— — любая полуобратная матрица к матрице S; 3) входные данные достаточно гладкие в U; 4) начальные и граничные данные удовлетворяют условиям A(x0,t)u(x0,t) = ^(t), B(x, t0)u(x, t0) = 0(x) и предполагаются согласованными в точке (0, 0) со своими производными. Тогда система уравнений (1) имеет решение в области U.
Доказательство. В силу леммы 2 после умножения системы (1) на матрицу P2 и замены переменной u = Q2z получаются три матричных уравнения с неизвестной вектор-функцией z = (zj zj zj)
z3 = £з,
где (gj gj gJ)T = Pf; Cj — блоки матрицы (7 = PADtQ + PBDxQ + PCQ. Из (15)
v(x,t + т) = v(x,t) + g(x,t) — т WAx + Lv(x,t) + h^K(x,jh,t)v(jh,t)
j=0
v(x,t0) = ^(x) — 0(t0), v(x0,t) = 0,
Dtzi + JD^zi + B13 D^z3 + (Ciizi + Ci2z2 = gi,
Dxz2 + C2izi + C22z2 = g2,
(15)
ж
жо
X _
где Ф(ж,£) = X(x,t)c(t) + J X(x,t)X-1(s, t)g2(s, t)ds; K(x, s,t) = —X(x,t)X-1(s,t) C21 (s, t),
0
X(x,t) — матрицант; c(t) — произвольная вектор-функция. Из (15) с помощью (16) получим начально-краевую задачу вида (14)
~ X
Dtz1 + JDXz1 + Си21 + J K(x, s, t)z1(s, t)ds = g1;
xo
21 (x, to) = ^1 (x), 21 (xo, t) = 01 (t) ,
где 71 = g1 — B13Dxg3 — C^(x,t); K(x,s,t) = C12K(x,s,t); P(x, 0^(x) = (^(x) (x))T ,
которая в силу леммы 3 имеет решение в области U. В (15) c(t) определяется однозначно: c(t) = X-1(x0,t)02(t) с помощью краевого условия z2(x0,t) = Ф2(t). Таким образом, показано, что решение существует. Теорема доказана. □
3. Анализ метода прямых
Одним из способов решения систем уравнений в частных производных является метод прямых. В исходной системе уравнений одну из частных производных заменяют конечной разностью, получая набор систем обыкновенных дифференциальных уравнений (ОДУ) [9].
В нашем случае вводим сетку по пространственной переменной xj = x0 + hj, j = 0,1, 2,..., N, h = (X — x0)/N и заменяем производную Dx конечной разностью (uj+1 — Uj)/h. При этом получаем набор задач Коши
Aj+1Uj+1 + Bj+1 [uj+1 — uj] /h + Cj+1Uj+1 = fj+1, (17)
где U0 = ^(t); Aj+1 = A(xj+1,t); Bj+1 = B(xJ+1,t); Cj+1 = С(xj+bt); fj+1 = f (xj+bt);
' = d/dt. Эту систему преобразуем к виду
Aj+1Uj+1 + Bj+1Uj+1 gj+1, u(xj+1, 0) 0(xj+1), (18)
где Bj+1 = B(xj+1,t)/h + С(xj+1, t); gj+1 = fj+1 + Bj+1Uj/h.
После дискретизации мы получили набор систем ОДУ с тождественно вырожденными матрицами при производной. Для решения этих систем нужны специализированные численные методы. Одним из таких методов является метод сплайн-коллокации [10]. Метод заключается в следующем. На каждом отрезке [tk , tk + тт] ищем приближение к решению x(t) в виде полинома m-й степени Lm(t), k = 0,1, 2,... (для одношагового процесса)
с неопределенными коэффициентами и требуем, чтобы этот полином удовлетворял в точках сетки решаемой системе. При этом мы получаем систему уравнений относительно коэффициентов полинома. Из соображений, связанных с обусловленностью, мы от системы относительно коэффициентов переходим к системе относительно значений полинома в точках сетки и таким образом получаем некоторую разностную схему.
Применим этот подход к системе
A(t)x(t) + B(t)x(t) = f (t), t G I = [t0,T]; (19)
x(t0) = x0, (20)
где A(t), B(t) — квадратные (n x п)-матрицы; det A(t) = 0 V t G I; f (t) — известная, x(t) — искомая n-мерные вектор-функции; x0 G Rn.
Соответствующая разностная схема (система линейных алгебраических уравнений СЛАУ) имеет вид
1 m
A(ifc + ¿т)- Yj,iXfcj + B(ifc + ¿т)xfc)j = f (tfc + ¿т), i = 1,m, (21)
"7” ' ^
T j=0
где к — номер шага процесса; (Т — ¿0)/тт < к < (Т — ¿0)/т, 7?,» — коэффициенты, отвечающие за аппроксимацию производных искомой вектор-функции ж(£). Они выбираются из условия
- m
p(ifc + ir)------y^7j,ix(tfc + ir )|| = O(rm).
rj=0
Формулы для вычисления 7j,i можно найти в [9]. Система (21) является замкнутой в том смысле, что число искомых векторов совпадает с числом уравнений. Решением системы (21) является блочный вектор (xfc,i, Xfc,2,Xfc,m), Xfc,0 считается известным. Поясним, что одной из особенностей метода является возможность на (к + 1)-м шаге численного процесса выбрать в качестве стартовой точки x0 любую компоненту блочного вектора
(xfc,1, xfc,2, •••, xfc,m) •
Сходимость схемы (21) докажем в предположении, что пучок матриц A A + B удовлетворяет критерию “ранг-степень” и в качестве стартовой точки берем Теорема 4. Пусть для задачи (19), (20) выполнены условия:
1) A(i), B(i), f (i) e Cm+1 (I);
2) rank A(i) = deg{det(AA(t) + B(i))} = const = d V i e I;
3) rankA(io) = rank{A(to), f (io) — B(io)xo}.
Тогда:
1) задача (19), (20) на отрезке I имеет единственное решение из класса Cm+1(/);
2) начиная с некоторого т < т*, СЛАУ (21) имеет единственное решение {xk>1, жк,2, xk,m}, для которого справедлива оценка ||x(tk) — xk;1|| = 0(тт).
Доказательство. В силу первого пункта утверждения теоремы задача (19), (20) имеет единственное решение x(i) e Cm+1(/), оно разложимо в ряд Тейлора с соответствующим количеством членов. Положим x(i) = Tm(i) + \P(i), ||Ф(*)|| = 0(тm), где Tm(i) — m-я частичная сумма ряда Тейлора, тогда, подставляя x(i) в уравнение (1), будем иметь
A(i)Tm(i) + B(i)Tm(i) = f(i) + e(i), 1Ш1 = 0(тm). (22)
С другой стороны, на каждом отрезке [ik,ik + тт] решение x(i) задачи (19), (20) аппроксимирует полином Lm(i) и выполняется равенство
A(tfc + ¿т )L m(ik + ¿т ) + B(ifc + ¿т )Lm(ik + ¿т ) = f (ifc + ¿т ). (23)
Находим разность между (22) и (23) и записываем разностную схему вида (21) относительно = ж (ifc + ¿т) — Xfc,j, где Xfc,j — значения многочлена из (23) в точках сетки.
Ф^ = т Ф — F, (24)
^ Afc,171,1 + тBfc,1 Afc,171,2 ... Afc,171,m 4
Afc,272,1 Afc,272,2 + т Bfc,2 ... Afc,272,m
Ф = . .
у Afc,m7m, 1 Afc ,m Tm,2 . . . Afc,m7m,m + тBfc,m у
( А,і7і,о N
272,0
^^k,0,
у Ак,т7т,0 у
где Ак,г, В^,г — матрицы, вычисленные в точке + ¿т. Умножим (24) слева на матрицу Р = diag{Pfcд, Рк,2,..., Рк,т} и сделаем замену переменной ш = 2у, 2 = diag{Qfcд, О^, 2к,т}.
Разложим диагональные блоки в матрице 2 в ряд Тейлора и, учитывая второе условие теоремы, а также равенство (10), перейдем к системам
Ук,% = 0, і = 1,т;
[Г + ** ] **+1 = — М* ** + Ф •
(25)
(26)
Здесь ||*к|| = 0(т), ||Ф|| = 0(гт+1), где векторы (25) и вектор гк+1 образуют вектор
Г = Е 0 7, 7
( 7і,і 7і,2
72,1 72,2
У 7т, 1 7т,2
71,т
72,т
7т,т
\
( *М \
Zk,2
\ *к,т у
М*
/" 7і,0Егі + *1 0 \
72,0Егі + *2 0
У 7т,0Егі + »т 0 У
||» II = 0(т)
(^' = 1,п, 0 — тензорное произведение матриц).
В силу свойств тензорного произведения и того, что матрица 7 невырожденная [10], матрица Г + ** обратима при достаточно малых Л Из (23) следует
**+1 = [Г 1 + ^*] ( —Мг*г + Ф)
где ||Ж*|| = О(Л). Поскольку мы берем в качестве начальной точки ж*д, запишем вычислительный процесс
**+1,1 = (Е 0 ••• 0) ^+1 • (27)
Нетрудно показать, что
/ Е + ^1 0 \
1 + $2 0
Г-1М = - . .
у + $т 0 У
где ||0г|| = 0(т)• Достаточно рассмотреть равенство
(7о 7Шт) = Т(т),
где Т(т) = (1 (1 + тт) • • • (1 + (тт)т) )Т , 7о = (71,0 72,о---7т,о)Т• Если умножить
„т1ТТТ„„ тт„ Л,-1
по-
следнее слева на 7 , то мы получим
Ііт 7 7о = -(1
т ^0
т
(28)
Значит, справедлива оценка ||^+1д|| < 5||^д|| + ||Д&||. По формуле геометрической прогрессии получаем оценку ||^д|| = 0(тт). Тогда по правилу треугольника, возвращаясь к старым переменным, запишем
Замечание 3. Из (28) следует, что в качестве стартового вектора можно выбирать не только ж*,1 для следующего шага, но и любой вектор из набора {ж*,2, ж*,3, •••, Жк,т}. Численные эксперименты показали, что выбор первой точки (особенно при применении метода сплайн-коллокации для нелинейных систем) обеспечивает большую численную устойчивость процесса вычислений.
Докажем следующее вспомогательное утверждение.
Лемма 4. Если в системе (1) пучок матриц АА + + С удовлетворяет двойному
критерию “ранг-степень” в области и, то в системе (18), начиная с некоторого к < к*, пучки матриц АА^+1 + В^+ при любом і удовлетворяют критерию “ранг-степень”.
Доказательство. Умножим пучок матриц АА^+1 + В^+ слева и справа на матрицы Р2 и ^2 соответственно и, учитывая лемму 2, получим
Выполняя элементарные преобразования над строками и столбцами выражения, стоящего в правой части, и учитывая, что при достаточно малых к < к* матрица £ /к + С22 обратима, получим
Следовательно, в силу леммы 1 пучок матриц АА^+1 + В.,+1 удовлетворяет критерию “ранг-
Теорема 5. Если для задачи (1), (2):
1) выполнены все условия теоремы 3;
2) на каждом слое х = ж0 + кі, і = 0,1, 2,...,Х, для задачи Коши (18) выполнены все условия теоремы 4;
3) начиная с некоторых т, к : т < т*, к < к*, шаги по пространственной и временной переменным связаны соотношением т = О (к),
то задача (1), (2) будет иметь решение щ*. и справедлива оценка ||м(ж^-,ік)— '•,*|| = О(к).
Доказательство. Докажем сначала сходимость метода прямых при предположении, что на каждом слое задача Коши решается точно. Положим г^+і = и*+1 — '-+1, где и*+1 = и*(ж,-+і,і), '+ — соответственно решения уравнений (1), (17). Подставляя в систему (17) решение системы (1) и вычитая из полученного соотношения саму систему (17), запишем
||х*,1 — < ||хМ — + ||х*,1 — = О(тт).
Теорема доказана.
□
степень.” Лемма доказана.
□
А7 + 1^7+1 + [В7+1/к + С7 + 1] Ч> + 1 = В7 + 1 /к"(,7 + Ь1'
(29)
где ||Ь1||с = 0(к). Покажем, что ||г-+1||с = О(к), ] = 0,1, 2,...,Ж Здесь и далее используется равномерная норма ||г-1| = шахЦг-,^, Ь £ [Ь0,Т], г = 1,п}, под нормами матриц понимаются нормы, согласованные с равномерной нормой.
В силу условия теоремы пучок матриц АА-+1 + ^В^+1/к + С-+1 удовлетворяет двойному критерию “ранг-степень”. Умножим уравнение (29) слева на матрицу Р2,^+1 и произведем замену переменной г-+1 = ^2,^+1^^+1, где Р2, ^2 — матрицы из леммы 2. Введем разбиение векторов и- = (г" г" г"-)Т , 61 = (V]1" V" . С учетом того, что = 0, получим
систему вида
£ 0\г +( Л'+1/к + ¿11 <12 ,
0 0) -+1 Су 21 £п-^/к + ¿22 У -+1
-+;/ £ ¿/л )+р2,-+1в-+1/к (о?- ^-+1) +р2,-+^^. (з°)
При достаточно малом к матрица £га-^/к+(722 обратима. Введем разбиение г- = (г""- г"-)Т. Умножим второе блочное уравнение в (30) слева на матрицу (£га_^/к + С22)-1 и, учитывая, что ||(£га_й/к + Су22)_1|| = О(к), перейдем к равенству
г2,-+1 = [£га-й + К2] г2,- + К1г1,- — К3г1,-+1 + 62, (31)
где К1; К2, К3 — матрицы, нормы которых имеют порядок О(к), норма вектора 62 порядка О(к2). Подставляя (31) в первое уравнение (30), получим линейное дифференциальное уравнение
¿і,і+і +
Л'+1/к + С11 + К4 ^1,і+1 = Зі(і) (32)
где Зі (і) = [7і+і/к + К5] гі;і +
Кз - С-
12
г2 - + 63, нормы матриц К, I = 4, 5, и вектора 63
порядка О (к). Выпишем решение уравнения (32)
£
г1-+1 = %(ф1-+1(Ьо) + J %(¿)%-1(з)д-(з)<^,
£о
где %(Ь) — матрицант уравнения (32) [11]. Напомним, что по условию теоремы 3 все собственные числа матриц /?+1 положительные, г-+1(Ь0) = 0. Из [12] известно, что
ехр ( J Атш(в)^ I < ||%(Ь)|| < ехр ( J Атах(«)^5
\£о / \£о
где Атщ(з) и Атах(з) — соответственно минимальное и максимальное собственные значения матрицы
Д(Ь) = — 2 [^1 + 1]) ^1 = ^-+1/к + + К4.
Теорема Гершгорина из [11] позволяет утверждать, что собственные числа матрицы Л(Ь) при достаточно малом к отрицательные и удовлетворяют оценке |Аг(Ь) — Аг/к| < с, с = ||Сп + К4||, Аг — собственные числа матрицы /-+1. Поскольку %(¿)%-1(з) является решением уравнения 21^ = — [/?+1/к + (711 + К4]г1;-+1 на отрезке [з,Ь], справедливо неравенство
||ё(¿)ё 1(в)|| < 1. Из уравнения Ё(¿) = (¿) получаем Ё 1 = — Ё 1Я1С учетом
последнего запишем
£ £
^ Ё(¿)Ё 1(в)з^- (з)^5 = —Ё(¿) ^ Ё 1Д]_ 1 (Д + Кб) £1,? + (Кз — ¿12)£2,.?' + Ьз ^5.
£о £о
Интегрируя равенство по частям и учитывая, что = 0(Л,), ||ё(¿)ё_1(^)|| < 1,
Ну Z(t)Z (s)C(s)ds|1 < C1||C||Ch,
¿0
где £(t) — некоторая дифференцируемая функция, получаем структуру переходного оператора для норм
|z1j+1||с ^ < ( 1+ c3h C4h ^ ( ||z1,j ||c І + b2
Z2,j+1 ||c J V C5h 1 + Cfefr y V ||Z2,j ||c
где ||&21| = O(h2). Отсюда следует, что ||zj+i||c < (1+ c7h)||z?||c + &2, где C7 = C3 + C4 + C5 + ce.
По формуле геометрической прогрессии имеем
||zj||с < [exp(c7(T - to)) - 1]62/c7h.
Так как ||b2||C = O(h2), то ||zj||C = O(h).
Обозначим через Uj^ решение уравнения (17), получаемое методом сплайн-коллокации, тогда в силу теоремы 4 и связи т = O(h ™ ) имеем
||u(xj,tfc) - uj-fc|| < ||u(xj,tfc) - Uj(tfc)|| + ||uj(tfc) - uj-fc|| = O(h) + O ((^) ) = O(h).
Теорема доказана. □
t
4. Результаты численных экспериментов
Для проведения численных экспериментов был разработан комплекс программ в среде Бе1рЫ 6.0. Он состоит из программы по решению систем обыкновенных дифференциальных уравнений методом сплайн-коллокации, двух программ по решению систем дифференциальных уравнений в частных производных методом прямых с дискретизацией по
переменным X и
В качестве тестовых примеров рассматривалось большое количество задач вида (1),
(2), удовлетворяющих условию (3), для которых выполнены все условия теоремы 3. Задачи подбирались так, что их точное решение было известно.
Пример 1.
/100\ / п 0 о \ / п п2 п3 \
А =10 0 0 , В = 0 п 1 ,С = 1 1 п2 ,
0 0 0 0 0 0 0 0 1
f = ( 0.1 (п + 2п2 + 2п3 + 3п4) 0.1(6п + 2п2 + Зп3) 0.3п )Т , U = [0,1] х [0,1], где п = ex+t.
Данная система удовлетворяет всем условиям теоремы 3 и имеет индекс (1,1). Известно ее точное решение и = (0.1^ 0.2^ 0.3^)Т . В процессе счета порядок сплайна задавался равным двум. Пространственный шаг Н варьировался от 0.1 до 0.0001, временной шаг выбирался соответственно условию теоремы 3. Получена абсолютная погрешность е = дН, где д « 0.8, е = тах ||и(х^, ) — и^||, ] = 1, 2,..., N к = 1, 2,..., N1.
Пример 2. Рассмотрим систему с коэффициентами
/1 0 0\ /А 0 0\ /а1 а2 а3^
А =|0 0 0] , В = 0 ех 1 , С = 0 7 0
\0 0 0/ \0 0 0/ \0 еж* 5
где А, 5, «¿, г = 1, 3, — числовые параметры, правой частью и краевыми условиями соответствующие решению и = еж+* (111) . Эта система гиперболическая и имеет индекс
(1, 2) при 5 = 1 и индекс (1,1) при 5 = 1. Элементарные вычисления показывают, что при применении метода прямых переходной оператор для второй компоненты имеет вид и2^+1 = ^(х, ¿)и2,^, где ^(х, ¿) = (е*х^+1 — )/(7Н), если 5 =1. Мы получаем неустойчи-
вый численный процесс, если ^(х,£)| > 1.
Были проведены эксперименты для систем индекса (1, 2), возникающих в некоторых физических приложениях. Система была получена заменой переменных из уравнения Кортвега —де Фриза — Бюргерса
I
2„ I „ п3 ^
Ар + «Ар + а2Лр + asB,p + ^ = -wJ Д,р
0
где а — некоторые константы. Она описывает распространение волн в парожидкостной среде [13]. Полученная система решалась методом прямых. Так как система является квазилинейной и содержит интегральный оператор, данные для расчета коэффициентов и интеграла брались с предыдущего пространственного слоя. Метод проверялся на системе с известным решением, а результаты расчета сравнивались с результатами физических экспериментов на высокотемпературном контуре Института систем энергетики им. Ме-лентьева СО РАН. В обоих случаях получены удовлетворительные результаты.
Список литературы
[1] Демиденко Г.В., УспЕнский С.В. Уравнения и системы, неразрешенные относительно старшей производной. Новосибирск: Науч. книга, 1998.
[2] Campbell S.L., Marzalek W. The index of infinite dimensional implicit system // Mathematical and Computer Modelling of System. 1999. Vol. 5, N 1. P. 18-42.
[3] Таиров Э.А., Злпов В.В. Интегральная модель нелинейной динамики парогенерирующего канала на основе аналитических решений // ВАНТ. Сер. Физика ядерных реакторов. 1991. Вып. 3. C. 14-20.
[4] Бояринцев Ю.Е., Чистяков В.Ф. Алгебро-дифференциальные системы. Новосибирск: Наука, 1998.
[5] Чистяков В.Ф., Щеглова А.А. Избранные главы теории алгебро-дифференциальных уравнений. Новосибирск: Наука, 2003.
[6] Чистяков В.Ф. Алгебро-дифференциальные операторы с конечномерным ядром. Новосибирск: Наука, 1996.
[7] Булатов М.В. Об одном семействе матричных троек // Ляпуновские чтения и презентация информационных технологий / ИДСТУ СО РАН, Иркутск, 2002. С. 10.
[8] Годунов С.К. Уравнения математической физики. М.: Наука, 1971.
[9] Березин И.С., Жидков Н.П. Методы вычислений. М.: Наука, 1966. Т. 2.
[10] Булатов М.В., Чистяков В.Ф. Применение коллокационных методов для решения сингулярных линейных систем ОДУ // Модели и методы исследования операций. Новосибирск: Наука. Сиб. отд-ние, 1988. С. 164-170.
[11] Гантмахер Ф.Р. Теория матриц. М.: Наука, 1967.
[12] Демидович Б.П. Лекции по математической теории устойчивости. М.: Наука, 1967.
[13] НАкоряков В.Е., Покусаев Б.Г., Шрейбер И.Р. Распространение волн в газо- и парожидкостных средах. Новосибирск, 1983.
Поступила в редакцию 17 мая 2004 г., в переработанном виде — 23 ноября 2004 г.
