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

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

Том 11, № 6, 2006

ОБ ОДНОМ ВАРИАНТЕ МЕТОДА ПРЯМЫХ ДЛЯ СИММЕТРИЧЕСКИХ t-ГИПЕРБОЛИЧЕСКИХ СИСТЕМ*

А. М. Блохин Институт математики СО РАН им. С. Л. Соболева,

Новосибирск, Россия e-mail: [email protected]

А. С. Ибрагимова Новосибирский государственный университет, Россия e-mail: [email protected]

A class of differential-difference models for symmetric systems is proposed. Approximate solutions for such models are studied.

Введение

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

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

1. Предварительные сведения

В качестве исходной математической модели возьмем смешанную задачу для симметрической t-гиперболической системы с диссипативными граничными условиями [1, 2].

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-00900 и № 06-08-00384) и междисциплинарного интеграционного проекта фундаментальных исследований СО РАН-2006 (грант № 46).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

Задача I. Найти решение и симметрической ¿-гиперболической системы (по Фридрих-су) с постоянными вещественными коэффициентами

A ■ Ut + B ■ Ux + C ■ Uy = 0 в области R+,

удовлетворяющее граничным условиям при x

0, l

U1 = U11

tII

S ■ U", x = R • U1

x

l

(1.1)

(1.2)

(1.3)

(1.4)

и начальным данным при t =0

U(0,x,y) = Uo(x,y), (x,y) G R,2.

Здесь

A = block diag(AI,AII,AI11);

B = block diag(/w0, —/w , On2 ) — блочно-диагональные матрицы; A1 = diag(ai,..., a^,);

aI1 = diag(aN0+1, . . . , aN0+Ni );

A111 = diag(aN0+N1+i,... , aw) — диагональные матрицы; a > 0, i = T^N, No + Ni + N2 = N;

/w0, Iw1 — единичные матрицы порядка N0, N1 соответственно; Ow2 — нулевая квадратная матрица порядка N2;

c

с1(= с1), с4(= с4), с6(= с6) —квадратные симметрические матрицы порядка N0, N1, N2 соответственно; с2, с3, с5 — матрицы размерности N0 х N1, N0 х N1 х N2 соответственно;

— вектор неизвестных функций,

ci c2 c

C2 c4 c

C3 c5 c

U = U(t, x, y) ( ui \

U1

( Uw0,1 \

U1

\uN0/

V0+1

U

iii

\uW0+Ni /

uN0+N1 + 1

\ Uw /

Л — вещественные постоянные матрицы размерности N0 х N1, N1 х N0 соответственно; Л+ = {(г,ж,у); ¿> 0, 0 < ж < /, у € Л1)}, Л2 = {(ж, у); 0 <ж</, у € Л1)}.

Замечание 1.1. Симметрическая система (1.1) с вещественными постоянными коэффициентами записана в так называемом каноническом виде (см. по этому поводу [1, 2]).

Будем полагать далее, что граничные условия (1.2), (1.3) являются строго диссипатив-ными [1], т.е.

— (B ■ U, U)|x=o = (U11, [/N1 — S* ■ S] ■ UII)|x=o > ko ■ (U11, U11) |x=o;

(1.5)

0

6

(В ■ и, и)|х=1 = (и1, [/N0 - я* ■ Я] ■ и1)|х=1 > к ■ (и1, и1) |л=г, (1.6)

где ко, кг > 0 — некоторые постоянные.

Умножая систему (1.1) скалярно на вектор 2 ■ И, можно на гладких решениях этой системы записать тождество интеграла энергии в дифференциальной форме

(А ■ И, И)4 + (В ■ И, И)х + (с ■ И, И)у = 0. (1.7)

Интегрируя затем тождество (1.7) по области Я2 и учитывая неравенства (1.5), (1.6), мы получим для Задачи I следующую априорную оценку (более подробно вывод априорной оценки описан в [1]):

3(Ь) < 3(0), Ь> 0. (1.8)

Здесь

3(Ь) = / /(А ^ И, И)^у,

Я2

3(0) = / /(А ■ Ио, Ио)^Ыу.

я?

При выводе оценки (1.8) мы полагали также, что

(И, И) ^ 0 при |у| ^ то.

Замечание 1.2. Теорема существования и единственности достаточно гладкого решения Задачи I со строго диссипативными граничными условиями может быть доказана с использованием априорной оценки (1.8) [1].

Рассмотрим так называемую дифференциально-разностную модель для нахождения приближенных решений Задачи I. При этом будем следовать принципу адекватности [3]: дифференциально-разностная модель для Задачи I будет конструироваться так, чтобы она допускала наличие дифференциально-разностного аналога априорной оценки (1.8). Наличие такого аналога априорной оценки будет означать также устойчивость предлагаемой дифференциально-разностной модели [3].

2. Дифференциально-разностная модель для Задачи I

Мы уже отмечали выше, что основная идея при построении дифференциально-разностной модели заключается в том, что в исходной задаче проводится дискретизация только по части независимых переменных (в частности, для Задачи I — по переменным Ь,у). Введем следующие обозначения:

И = И(к ■ Ну ) = И* (х) = И,- (х), 0 < х < /, к, Ц | = 0,1,...,

И=ип*>. <=

9 - 9-1 г 9 + 9-1

По = , ьу = ~2~

— разностные операторы, 9, 9-1 —операторы сдвига: 9±1И^-(х) = И^±1(ж), 9+1 = 9, Д, Ну — шаги разностной сетки по переменным ¿, у.

Сопоставим Задаче I такую дифференциально-разностную модель.

Задача II. Ищем вектор И^(ж), удовлетворяющий следующим соотношениям:

А • ИИ = А • Ьу И -ДВ • (И - г у с • П0И; (2.1)

0 < ж < /, к, |7| = 0,1,...; И = 5 • И?, ж = 0, к, 7| = 0,1,...; (2.2)

И1 = Л • И;, ж = /, к, 7| = 0,1,...; (2.3)

И0(ж) = И0(ж,7 • йу), 0 < ж < /, Ц1 = 0,1,... (2.4)

Здесь гу = .

Замечание 2.1. Равенство (2.1) получается так. В системе (1.1) аппроксимируем производные И Иу такими разностными выражениями, как

И;(ж) - И;(ж) П0И (ж)

д , ^ .

В результате система (1.1) переписывается в виде соотношения (2.1).

С помощью техники построения дифференциально-разностного аналога диссипатив-ного интеграла энергии для Задачи I [3] докажем корректность (устойчивость) Задачи II. Для этого перепишем систему (2.1) так:

А • "И = • И;+1 + ^2 • И-1 - ДВ • С ИИ, (2.1')

где

Di = 1(A - ryc), D2 = 2(A + ryc).

Систему (2.1') умножим скалярно на вектор и:

(А • ИИ, ИИ) = р! • И;+1, ИИ) + (а • И;-1, ИИ) - Д(В • с ИИ, ИИ). (2.5)

Заметим, что матрицы А12 > 0, если

< ^0.

Здесь v0 = min а^; ||с|| = max А — операторная норма матрицы c; Адс), i = 1,N, i=1,N i=1,N

собственные числа матрицы с. Поскольку

(Di ■ Uj+i, U) < 2(Di ■ Uj+i, Uj+i) + 1(Di ■ U, U),

(D2 ■ Uj-i,U) < 1(D2 ■ Uj-i,Uj-i) + 1(D2 ■ U,U), из (2.5) следует неравенство

(A ■ U, U) + АС(B ■ U, U) < (Di ■ Uj+i, Uj+i) + (D2 ■ Uj-i, Uj-i). (2.6)

Обе части неравенства (2.6) умножим на Ну, просуммируем полученное выражение по Щ от —то до то, затем проинтегрируем его по х от нуля до /:

Зк+1 + ДЛу £ {(В ■ и, и)|х=1 — (В ■ и, и)|*=о}< Зк, (2.7)

j=-x

где

I

X

Зк = Ну ^ / (А ■ ик (х), ик (х))йх.

0

Наконец, учитывая строгую диссипативность граничных условий (см. неравенства (1.5) и (1.6)), из неравенства (2.7) получаем дифференциально-разностный аналог априорной оценки (1.8):

Зк < Зо, к = 1, 2,..., (2.8)

означающий устойчивость дифференциально-разностной модели в энергетической норме

\Пк [3].

Замечание 2.2. При выводе неравенств (2.6)-(2.8) мы воспользовались также почти очевидными соотношениями

(В ■ с и, и ) = 1 С (В ■ и, X}),

(

{(А ■ и^+1, и^+1) + (Д2 ■ и-1, и^-1)} = £ (А ■ и, и).

{(^1 ■ и+Ь + (^2 ■ }j-1, }j-1)

j = — X j= — (X

Замечание 2.3. Без труда можно предложить и другие варианты дифференциально-разностных моделей в случае Задачи I [3], для которых можно получить дифференциально-разностный аналог априорной оценки (1.8).

3. Разрешимость Задачи II

Рассмотрим теперь вопрос о нахождении вектора и(х) на (к + 1)-м слое по известному вектору ик(х) на к-м слое. С учетом специфики матриц А, В (см. разд. 1) систему (2.1) можно переписать так:

сV = А V + Т, (3.1)

и111 = Ьуи111 — Гу(Аш)-1сш ■ пои, (3.2) 0 < х < /, к, |Щ| = 0,1,...

и1 и1

Здесь V = Ы; Ы = ( тт11 I ,

А = 1 Ь1оск diag(—А1, А11), А = Ь1оск diag(AI, А11), В = Ь1оск diag(/N0, —) — блочно-диагональные матрицы;

Т = 1 В(АЬу Ы — Гу с пои),

с111 = (с с5 с6) — блочная матрица размерности N х N; с = ( с с ^ ) — блочная

УС2 С4 С5)

матрица размерности (N0 + х N.

Граничные условия для систем обыкновенных дифференциальных уравнений (3.1) примут следующий вид (см. соотношения (2.2), (2.3)):

(/N0, -5) • У(0) = 0; (3.3)

(-Л, /^) • У(/) = 0. (3.4)

Здесь (/^, -$), (-Л, /^) — блочные матрицы размерности N0 х (N0 + и N х (N0 + N1) соответственно.

Поскольку (3.1) — это система обыкновенных дифференциальных уравнений с постоянными коэффициентами, легко можно найти решение краевой задачи (3.1), (3.3), (3.4) (см. [4, 5]). Ниже мы изложим технологию нахождения приближенных решений этой краевой задачи, основанную на технике сплайн-функций. Заметим, что этот способ пригоден и для случая переменных коэффициентов исходной системы (1.1). Будем искать решение краевой задачи (3.1), (3.3), (3.4) приближенно в виде кубического нелокального сплайна класса С1 (см. по этому поводу [6]):

Яс(ж) = (1 - 0)2(1 + 2^(гйж) + 02(3 - 2^((г + 1)й*) +

+М(1 - 0)[ш,(1 - 0) - шт0], (3.5)

x е [ihx, (i + 1)hx], i = 0, J — 1, hx J = /.

Здесь $ = (x — ihx)/hx; m, = (Uj|ж=а*•

Кубический сплайн (3.5) непрерывен вместе со своей первой производной всюду на отрезке [0, /]. Следуя [6], будем полагать, что вторая и третья производные от кубического сплайна (3.5) удовлетворяют соотношениям

S" (ihx + 0) — S" (ihx — 0) =

а, ■ {S£(ihx + 0) — S'("/(ihx — 0)}, i = 1, J — 1, (3.6)

где а, — заданные числа.

Вычисляя с помощью (3.5) векторы

s"(ihx + 0), S"(ihx — 0), S'("/(ihx + 0), S"^ — 0)

и подставляя их в (3.6), получаем

£

(1 — 3а,) ■ m,-i + 4 ■ m, + (1 + 3а,) ■ m,+ = 3(1 + 2а,)^Uj (ihx) +

+3(1 - 2а,)(гйл), г =1,7 - 1. (3.7)

Здесь а, = а,/Л,ж; £ = ("0 - 1), £ = (1 - "0-1) — разностные операторы; — операторы

сдвига: (гЛ,ж) = Ц0((г ± 1)Л.Ж).

Рассматривая равенство (3.7) на (к+1)-м слое и исключая вектор V с помощью системы (3.1), мы придем к такому выражению:

{(1 - 3а,)А- 3(1 - 25,)гжВ} • М,-1 + (4А - 12а,ГхВ) • М,+

+{(1 + 35г)А + 3(1 + 25г)гжВ} ■ Мт =

£ £ _

= 3(1 + 25,)^Z + 3(1 - 2а,)7, г = 1,3 - 1, (3.8)

' ¿х 'х

где Гх = А/'ж; 7 = - Гус ■ П0И М, = (V|ж=Ах•

Для нахождения производных М,, г = 0, 3, с помощью системы алгебраических уравнений (3.8) необходимо задать краевые условия при г = 0 и г = 3 .С этой целью положим в системе (3.1) х = 0 и х = /. В итоге получим

Мо = А( ^ ) ■ И (0)+ Т(0); (3.9)

М^ = А^■ и1 (/)+ Т(/). (3.10)

Здесь ( ) , ( ^ ) — блочные матрицы размерности (N0 + Ж1) х N и (N0 + Ж1) х N

\ Я , соответственно.

Для окончательного определения агрегатов М0, MJ к (3.9), (3.10) добавим условия вида [1, 3]

Ии11 (0) = ип('х); (3.11)

х}1(/) = и!(/ - 'х). (3.12)

Система (3.8)—(3.10) представляет собой замкнутую систему алгебраических уравнений относительно неизвестных М,, г = 0, 3 (с учетом условий (3.11), (3.12)). Решение ее можно найти, например, методом прогонки [7, 8]. После этого значения решения V(гhx), г = 1,3 - 1, на (к + 1)-м слое находятся с помощью системы (3.1). Параметры а,, г =1,3 - 1, гх выбираются из условий устойчивости прогонки.

В нашем случае эти условия достаточно грубо можно описать так. Обозначим п-ю компоненту агрегата М,, г = 0,3, 1 < п < N + N1, через и,. Тогда, учитывая структуру матриц А, В, из (3.8)—(3.10) получаем разностную краевую задачу

_ и0 = ^0,

свд— + + с,и,+1 = г =1,3 - 1, } (3.13)

^ = ^.

Здесь

а, = (1 - 3а,)а„ - 3(1 - 2а,)гж^„;

с, = (1 + 3а,)а„ + 3(1 + 2а,)гх^га; = ±1.

Агрегаты г = 0,3, легко могут быть выписаны. Выбирая, например, числа а,, г = 1,3 - 1, так, чтобы

1 „ 1

- 3 < а, < 3, г =1,3 - 1, (3.14)

а шаг А был достаточно малым, мы придем к выводу, что выполнены условия хорошей обусловленности (устойчивости) разностной краевой задачи (3.13) (см. [7]):

|С| > |с| + |с,|, г = 1,3 - 1.

4. Сведение алгоритма с использованием техники сплайн-функций к разностной схеме

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

и = И(кД, г'х, ^) = ик = И,, = и, = и,-

— сеточная вектор-функция, к, ^| =0,1,... , г = 0, 7. Система (3.1) при ж = гЛ,х запишется так:

М, = А ■ V, + г = 1,7 - 1. (4.1)

Здесь

= -АЬуЦ, - Вс ■ ——- = — -В ■ 7,, 'у Д

7, = АЬуЦ, - г у с ■ П0И,. Подставляя вектор М, из (4.1) в (3.8), получим

{(1 - 3а,)А - 3(1 - 2а,)гхВ}А ■ V,-l + (4А - 12а,ГхВ)А ■ V,+ +{(1 + 3а,)А + 3(1 + 2а,)гх В}А ■ v,+l + {(1 - 3а,)А - 3(1 - 2а, )гжВ} ^ В ■ 7,-1+

+ (4А - 12а,гжВ4)^В4 ■ 7, + {(1 + 3а,)А + 3(1 + 2а,)гхВ}^В4 ■ 7,+1 =

и , о,-,

= 3(1 + 2а,)7, + 3(1 - 2а,)7,, г = 1,7 - 1. (4.2)

' X ' X

С учетом структуры матриц А, В, А от (4.2) после несложных, но громоздких выкладок приходим к такому выражению:

Ат^и,. + (1 + а,) гжВ ■ + (1 - ал ^В ■ + гуС ■ пь^И,- = 0,

к, Ц | = 0,1,..., г =1,7 - 1, (4.3)

где

, (1 - 3а,)^-1 + 4 + (1 + 3а,)^ т = (<£ - Ьу), ^ =-6-

— разностные операторы; ^ — оператор сдвига: = Ик.+1 = "И,,; = V,,; остальные

обозначения даны выше. Наконец, систему (3.2) перепишем в такой форме:

А111 ■ т^И, + гус111 ■ П0^И , = 0, к, Ц| = 0,1,..., г = 1,7 - 1. (4.4)

Объединяя системы (4.3), (4.4), получим (по терминалогии из монографии [3]) явно-неявную разностную схему для смешанной Задачи I:

А ■ т^хИк,- + ( 1 + ал гхВ ■ ,, + (1 - ал гхВ ■ £И,,- + гу с ■ пь^хИк,- = 0, (4.5)

к, Ц | = 0,1,..., г = 1,3 - 1.

Для определения сеточной вектор-функции и мы добавим к разностной схеме (4.5) граничные условия при г = 0 и г = 3 (см. условия (1.2), (1.3), (2.2), (2.3) и условия (3.11), (3.12)):

и 0; и Ъ и

S • И ,

я • и

и11

ъ ,

И Ъ = -1,;,

к, Ц| = 0,1,...

Векторы "И"ИЪ; определяются с помощью системы (3.2). При к чальные данные (см. условия (1.4) и (2.4)):

(4.6)

0 выполняются на-

и:

И0(г^ж, ^'Л-у), г = 0, 3,

0, 1,

(4.7)

Решение разностной смешанной задачи (4.5)-(4.7) находится методом прогонки. Относительно чисел с/г, г =1,3 — 1, будем полагать, что выполнены неравенства (3.14).

Замечание 4.1. Используя технику интегралов энергии, разработанную в [3], рассмотрим вопрос об устойчивости разностной схемы (4.5).

Замечание 4.2. Можно предложить и другие варианты реализации граничных условий (вместо (4.6)) для разностной схемы (4.5) (см. по этому поводу [3] и разд. 5).

о

5. Устойчивость явно-неявной разностной схемы (4.5)

Перепишем разностную схему (4.5) так:

А • ^ + 2г* • В • = А • ^"и^1 + • (5.1)

Здесь = {(1 + 25г)£ + (1 — 2с/г)£} — разностный оператор; матрицы Д12 описаны выше.

Следуя рекомендациям из [3], в качестве граничных условий возьмем такие соотношения:

И0 = • И0 — 7 • И1,

И01 = и?,

И = и -1,

(5.2)

и = АЯ • иъ — 71 • иЪ-1,

к, Ц | = 0,1,...

Здесь в, въ 7, 71 > 0 — вещественные параметры, причем 1 + 7 = в, 1 + 71 = А.

Умножим систему (5.1) скалярно на затем — на и просуммируем полученное выражение по г от 1 до 3 — 1, по ] — от —то до то (см. разд. 2). В итоге запишем

Ъ —1 х

Л+1 + £ (В • , V*и) < 3к, (5.3)

г=1 ;=—х

где

Ъ — 1 х

3*; = ^ ^ (А • ^и, ^и).

г=1 ;=—х

Рассмотрим случай, когда все а, = 0. Тогда

^х = £ + £, V* = ^ + + 4). 6

Обозначим через К следующий агрегат:

1 J-1 1 J-1

К = 6 Е(В(£ + ё)И, (^ + Г1 + 4) И) = - - Г!)И, (^ + Г1 + 4)И) =

,=1 ,=1

1 J-1 ^^ ^ ^ Л Л __

^{(В ■ ^И) - (В ■ ^-1И,^-1И) +4(В ■ и) - 4(В ■ ^-1И, и)} =

- ■ 1

,=1

1 {-(В ■ И0,ИИ0) - (В ■ ИИ 1,ИИ 1) + (В ■ ИJ-1,ИИJ-l) + (В ■ ИИJ,ИJ)-6

-4(В ■ И0, ИИ 1) + 4(В ■ ИИJ-1, ИJ)} = 1 {К1 + К2}.

6

Здесь

К1 = -(В ■ и0, и0) - (В ■ и 1, и 1) - 4(В ■ и0, и 1), К2 = (в ■ ИИ J-1, И J-l) + (В ■ И J, И J) + 4(В ■ И J-l, ИИ J).

Поскольку в силу специфики матрицы В

К1 = -(И0,И0) + (И01,И01) - (И 1,И1) + (И 11,И11) - 4(И0,И1) + 4(И01,И11),

учитывая граничные условия (5.2), получаем

К1 = 6(И 11, И11) - (И 1, И1) - [72 ■ (И 1, И1) - 2вт ■ (И 1, 5 ■ И11) + в2(5 ■ И11, 5 ■ И11)]-

-4[в ■ (И 1,5 ■ И11) - 7 ■ (И 1, И1)] =

= 6(И 11, И11) + (-1 + 47 - 72) ■ (И 1, И1) - в2(5 ■ И11, 5 ■ И11) - 2в(2 - 7) ■ (И 1,5 ■ И11). (5.4)

Из неравенства (1.5) следует

> (к0 - 1)/м1. (5.5)

Кроме того, справедливо очевидное неравенство (неравенство Шварца — Коши [2]):

2|(И 1,5 ■ И 11)| < (И 1, И 1) + (И^5 5 ■ ИИ 11). (5.6)

Учитывая (5.5), (5.6) и полагая для определенности -1 < 7 < 2, из (5.4) получаем

К1 > 40 ■ (И 11, ИИ 11) + 40 ■ (И 1, И1), (5.7)

где

40 = 3^0 = 3[(1+ 7)к0 - (7 - 1)], 4 = 3^0 = 3(7 - 1),

7 - 1

причем д0 > 0, д0 > 0, если 1 < 7 < 2, к0 >

Аналогично для агрегата К2 выводим неравенство

К > / • (иЪ—1, х}Ъ—1) + / • (хиЪ—1, тиЪ—1), (5.8)

где

/ = 3® = 3[(1 + 71)кг — (71 — 1)], & = 30г = 3(71 — 1),

71 — 1

причем ф > 0, о > 0, если 1 < 71 < 2, кг > -.

71 + 1

С учетом неравенств (5.7), (5.8) неравенство (5.3) примет вид

3к+1 < 3к, к = 0,1,...

Здесь

J = Jk + 2Ah^ {qo ■ (U?, U?) + go ■ (Ul, U?) + qi ■ (UJUj-i) + gi ■ Uj-i)}.

Значит, разностная модель (5.1), (5.2), (4.7) устойчива в "энергетической норме" \/3к.

Замечание 5.1. Мы взяли слова "энергетическая норма" в кавычки, поскольку, строго говоря, еще требуется доказать, что

J > 0 при Uj = 0.

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

[1] Годунов С.К. Уравнения математической физики. М.: Наука, 1979.

[2] Блохин А.М. Элементы теории гиперболических систем и уравнений: Учеб. пособие. Новосибирск: Изд-во НГУ, 1995.

[3] Блохин А.М., Алаев Р.Д. Интегралы энергии и их приложения к исследованию устойчивости разностных схем. Новосибирск: Изд-во НГУ, 1993.

[4] Годунов С.К. Матричная экспонента, матрица Грина и условие Лопатинского: Учеб. пособие. Новосибирск: Изд-во НГУ, 1983.

[5] Белман Р. Введение в теорию матриц. М.: Наука, 1976.

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

[7] Годунов С.К., Рябенький В.С. Разностные схемы. М.: Наука, 1973.

[8] Бабенко К.И. Основы численного анализа. М.; Ижевск: НИЦ "Регулярная и хаотическая динамика", 2002.

Поступила в редакцию 12 мая 2006 г.