CC BY
18ВЕСТН. МОСК. УН-ТА. СЕР.1, МАТЕМАТИКА. МЕХАНИКА. 2009. №5
11
УДК 519.6, 517.9
ДВАЖДЫ НЕПРЕРЫВНО ДИФФЕРЕНЦИРУЕМЫЙ ПОЛУЛОКАЛЬНЫЙ
СГЛАЖИВАЮЩИЙ СПЛАЙН
Д. А. Силаев1
Рассматриваются дважды непрерывно дифференцируемые периодический и непериодический полулокальные сглаживающие сплайны, или ¿"-сплайны класса C2, состоящие из полиномов пятой степени. Первые три коэффициента каждого полинома задаются значением предыдущего полинома и его двух производных в точке склейки, коэффициенты же при старших степенях полинома определяются по методу наименьших квадратов. Эти условия дополнены условием периодичности сплайн-функций на всем отрезке определения или начальными условиями, доказаны теоремы существования и единственности, установлены условия устойчивости и сходимости таких сплайнов.
Ключевые слова: аппроксимация, сплайн, сглаживание, полулокальность, полином, пятая степень, численные методы, математическая физика.
Twice continuously differentiable periodic local and semilocal smoothing splines, or ¿-splines from the class C2 are considered. These splines consist of polynomials of 5th degree, first three coefficients of each polynomial are determined by values of the previous polynomial and two its derivatives at the point of splice, coefficients at higher terms of the polynomial are determined by the least squares method. These conditions are supplemented by the periodicity condition for the spline function on the whole segment of definition or by initial conditions. Uniqueness and existence theorems are proved. Stability and convergence conditions for these splines are established.
Key words: approximation, spline, smoothing, semilocality, polynomial, fifth degree, numerical methods, mathematical physics.
Введение. Данная работа посвящена традиционной задаче восстановления значений функции и производных этой функции в промежуточных точках по заданной таблице. Эта задача эффективно решается с помощью построения тем или иным способом кусочно-полиномиальной функции, или сплайн-функции, приближающей исходную функцию по ее таблице значений в заданных точках. Полулокальные сглаживающие сплайны, или ¿-сплайны, были введены Д. А. Силаевым [1—5]. В отличие от ранее рассмотренных непериодических S-сплайнов класса C2 [2] в настоящей работе изменена схема построения таких сплайнов, что привело к существенному облегчению условий их устойчивости. Устойчивость периодического сплайна, в частности, проявляется в "притяжении" непериодических S-сплайнов, отличающихся различными начальными условиями, к периодическому ¿-сплайну. Указанный метод обеспечивает сглаживание результатов измерений и может быть использован при автоматизированной обработке информации.
1. Построение системы уравнений для S-сплайна. Рассмотрим на отрезке [a,b] равномерную сетку Xk = a + kh, k = 0,...,K, h = (b — a)/K — шаг сетки. Разобьем отрезок [a, b] на группы, для чего введем еще одну равномерную сетку = a + lH, H = mh, m E N. Таким образом, переходя из одной группы в другую, мы осуществляем сдвиг системы координат и рассматриваем каждый l-й полином на отрезке [0,Н]. Пусть y = (yo,yi,...,yk) E R^+1 — значения приближаемой периодической функции на этой сетке. Обозначим через
P'S = I u : u(x) = a0 + a1x + a2x2 + ^ aiX% \
^ i=3 J
множество полиномов степени n с фиксированными коэффициентами ao,ai,a2. Рассмотрим функционал
1 Силаев Дмитрий Алексеевич — канд. физ.-мат. наук, доцент каф. общих проблем управления мех.-мат. ф-та МГУ, e-mail: dasilaev@mail.ru.
12
ВЕСТН. моек. УН-ТА. СЕР.1, МАТЕМАТИКА. МЕХАНИКА. 2009. №5
м
ф1(и) + кН) - Ут1+к)2
к=0
В классе РЩ ищется такой полином до, который минимизирует функционал
м
ф1(и) = + кН) — ут1+к)2 — ш1п(аэ, аА ,...,ап)
„ ■ к -^
к=0
и удовлетворяет следующим условиям:
а10 =д1-1(£1 -6-1) =91-\{Д), а[=д'1_1(Н), а12 = ^д'1_1{Н) при 1 = 1,...,Ь-1. (1)
Так как а) = до(0),а1 = д1 (0),а2 = д1 (0)/2, то условия (1) есть условия гладкой склейки двух последовательных полиномов. В непериодическом случае начальные коэффициенты а0,а1,а2 задаются начальными условиями Уо,Уо,Уо/2?. В периодическом случае условие (1) при I = 0 есть условие периодичности 5-сплайна (считаем д-\ = дь-1 ). Здесь Ь — число групп, на которые разбита исходная таблица значений приближаемой функции, или число полиномов, составляющих сплайн. Кроме того, здесь М + 1 — количество точек осреднения; т + 1 — количество точек, входящих в область определения 1-го полинома до; — точка привязки полинома до; М — т + 1 — число таких точек, значения которых участвуют при определении двух соседних полиномов, составляющих 5-сплайн, М ^ т + 1. Можно предполагать, что значения заданной функции Ук известны с некоторой точностью, например они есть результаты каких-либо измерений. Будем предполагать тогда, что с уменьшением шага Н будет увеличиваться точность измерения, а именно будем предполагать, что если периодическая функция / £ С6 [а, Ь] задана в узлах равномерной сетки Хк = а + кН,к = 0,...,К, своими значениями yk , то \ук — f(хк)\ ^ СН6+£,е ^ 0. В дальнейшем степень полинома п = 5.
Определение. 5-сплайном назовем функцию вЩ, м(х), которая совпадает с полиномом до(х) на отрезке ^ ^ х < &+1.
Будем минимизировать функционал Фо по коэффициентам аз ,а4 ,а5. Для этого Фо(д) продифференцируем по этим коэффициентам и приравняем к нулю. Получим
м м мм
.3^7,6 , „о 1,4^7,7 , А _ . „о „¡и „о г,21„2\,„3П
к6 + а14Н4^ к7 + а15Н5 ^ к8 = ^[(Уто+к — а10 — а\Нк — а12Н2к2)к3],
к=0 к=0 к=0 к=0
м м мм
к1 + а14Н4^ к8 + а15Н5 ^ к9 = ^[(Уто+к — а) — а\Нк — а12Н2к2)к4], (2)
к=0 к=0 к=0 к=0
м м мм
к8 + а14Н4^ к9 + а15Н5 ^ к10 = ^[(Уш1+к — а) — а\Нк — а12Н2к2)к5].
к=0 к=0 к=0 к=0 Введем обозначения:
мм
~
= ^ к3, 3 = ^}(Ушо+к — а) — а\Нк — а^Н2к2)к3+2], (3)
к=0 к=0
учитывая которые, перепишем систему уравнений (2) в виде
а13 Н3в6 + а14Н4Б7 + а15Н5Б8 = е\, а3 Н3Б1 + а14Н4Б8 + а15Н5Б9 = е\, (4)
. а3Н3в8 + а14Н4в9 + а15Н5в10 = е33.
2В случае, если функция задана таблицей, у0,у0 можно определить с помощью формул численного дифференцирования
высокого порядка аппроксимации, например уо = — ^(147г/о — 360г/1 + 450г/2 — 400?/з + 225г/4 — 72уъ + 10уь)/Н + о(/гБ), уо =
1 180
(812уо - 3132у1 + 5265у2 - 5080уз + 2970у4 - 972у5 + 137уа)/к2 + о(к4).
Система линейных алгебраических уравнений, которой должны удовлетворять коэффициенты полиномов ¿-сплайна, состоит из уравнений двух видов: (а) уравнения склейки для каждой пары последовательных полиномов (1); (б) уравнений для определения коэффициентов при старших степенях полиномов по коэффициентам при нулевой, первой и второй степенях. Сделаем замену с^ = сцНг, г = 0,1,..., 5. При этом уравнения (а) примут вид
а.
I-1
1
,2т, 1-1
1-1
1-1
,5л1-1
0 + та{~ + т2 а2_ + т3а3_ + ш4а4~ + ш5а5~ а1-1 + 2т а,1-1 + 3т2 С1-1 + 4т3 а4-1 + 5т4 а5~1 = а\, а- + 3т а3-1 + 6т2 а4-1 + 10т3 а15-1 = а\.
(5)
Подставив выражения с) из (3) в систему (4), получим уравнения (б):
¿3 а0 + ¿4 а1 + ¿5 а2 + ¿6 а3 + ¿7 а4 + ¿8 а5 = Р\, ¿4 а0 + ¿5 а1 + ¿6 а2 + ¿7 а3 + ¿8 а4 + ¿9 а5 = р2, ¿5 а0 + ¿6 а1 + ¿7 а2 + ¿8 а3 + ¿9 а4 + ¿10 а5 = р3 ,
где
м
Р) = Е Ут+к к3+2.
к=0
Здесь I = 0,...,Ь — 1 — номер полинома, причем если I = 0, то в периодическом случае а1 1 означает
ггЬ-1
. Перенося в уравнениях (5) С0, а!, а2 в левую часть, получим систему уравнений для определения всех коэффициентов полиномов, составляющих периодический ¿-сплайн. В дальнейшем, если это не вызовет путаницы, мы будем опускать волну над переменными к.
Запишем полученную систему в матричной форме. Для этого обозначим:
¿3 ¿4 ¿5 ¿6 ¿7 ¿8 1 т 2 т2 т3 т4 5 т5
А1 = ¿4 ¿5 ¿6 , ¿2 = ¿7 ¿8 ¿9 , В1 = 0 1 2т , В2 = 3т2 4т3 5т4
¿5 ¿6 ¿7 ¿8 ¿9 ¿10 0 0 1 3т 6т2 10т3
Пусть, кроме того,
р
X2 =
X 21+1 =
где I = 0,1,. ..,Ь—1. Тогда систему уравнений для определения коэффициентов периодического ¿-сплайна можно записать в виде
/-Е 0 0 0
¿1 ¿2 0 0
В1 В2 —Е 0
0 0 А1 ¿2
0 0 В1 В2
V 0 0 0 0
В1 0 0 0 0
В2 0 0 0 0
А1 А2/
( X0 \
X1 X 2 X 3 X 4
XX 2Ь-1)
0
р0
0 р1
0
\р ь-1)
(6)
1 0 0 0 0 0
Матрицу этой системы обозначим О. Здесь Е = 0 1 0 ,0= 0 0 0
0 0 1 0 0 0
единичная и нулевая матрицы.
Отличие от аналогичной системы для непериодического ¿-сплайна, рассмотренного в работе [2], заключается в замене первой строки в (6) на стартовые условия, которые можно записать следующим образом:
(Е 0 0
имеет вид
Ш =
У°,
0)Х° = У°, где У° = | Ну° | . В работе [2] показано, что определитель матрицы
898128000МЗ ^М ~ 2^ + + 3) (5ММ + З5м13 + 25М12 - 305М11 - 630М10 + 480М9 + 2050М8+
+535М7 - 3010М6 - 4030М5 + 2370М4 + 8685М3 + 1278М2 - 3312М + 864)(М + 2)2(М - 1)2(М + 1)3.
Полученный полином имеет три положительных вещественных корня {1, 56081; 2; 2, 55841}, поэтому при М ^ 3 существует Л--1.
Система (6) распадается. Преобразуем нечетные строки матрицы С так, чтобы исчезла зависимость в (6) от нечетных неизвестных X 1,Х3,..., Х2Ь-1. Для этого из третьей строки вычтем вторую, умноженную на матрицу В2(Л2)-1. Тогда первые три строки будут иметь вид
-Е 0 0 0 Л1 Л2 0 0 В1 - В2Л-1 Л1 0 -Е 0
0
р °
-В2Л-1Р °
Аналогичные преобразования выполним для каждой пары соседних строк с номерами 21 + 1 и 21, где I = 0,1,... ,Ь - 1, а также для первой и последней строки. Обратим внимание на то, что получающаяся в результате таких преобразований матрица будет иметь вид и = В1 - В2Л-1 Л1. После преобразований система (6) расщепляется на систему для отыскания X2, I = 0,1,... ,Ь - 1, вида
(7)
(- -Е 0 0 . .. и \ / X ° \ /- -В2Л-1 р ь-1\
и -Е 0. . 0 X 2 -В2Л-1Р °
0и -Е. . 0 X 4 = -В2Л-1Р1
V 0 0 0 . ..и- Е ^2Ь- 2 V -В2Л-1 рь-2)
и следующие уравнения для отыскания X21+1:
Л2Х21+1 = Р - Л1Х21, I = 0,1,...,Ь - 1.
Найдем выражение для матрицы и. Введем обозначения:
(8)
Бг Б, Б к ^
= ^ | Бг+1 Бк+1 | , = ёе1(Л2).
уБг+2 Б,+2 Бк+2/
Заметим, что ¿е!(Л2) = Тб78. Тогда
Л-1 Л1 =
и элементы матрицы и будут иметь вид
пп = 1 - гз78т3 - гб38т4 - 1673т5, П1з = т2 - г578т3 - гб58т4 - 1б75>тъ, П22 = 1 - 31478т2 - 4Ьб48т3 - 5Ьб74т4, и31 = -31378т - 61638т2 - 10£673т3,
tз78 t478 t578 t638 t648 t658 , t673 ^674 ^675 ,
и12 = т - 1478т3 - 1648т4 - 1674т5, П21 = -3г378т2 - 41638т3 - 51673т4, и23 = 2т - 31578т2 - 4Ь658т3 - 5Ь675т4, и32 = -31478т - 61648т2 - 10^74т3,
и33 = 1 - 3г578 т - 6г658 т2 - 10£675 т3.
Легко заметить, что и^ = иЪ] = ^¿^"(^.у), 3 = 1)2,3. Отличие от формул (8) работы [2] состоит
в том, что здесь выписана система для определения коэффициентов полиномов, составляющих сплайн, а в работе [2] в качестве переменных использовались дг (0) = а0, дг (0) = а1, д1 (0) = 2а^.
2. Существование и единственность ¿-сплайна. Теорема 1. Пусть числа т и М ^ 3 таковы, что собственные числа матрицы и не равны корню степени Ь из единицы (здесь Ь — число полиномов, составляющих сплайн). Тогда для любой периодической функции /(х), заданной на отрезке [а,Ь] своими значениями Ук в точках Хк = а + кН, Н = (Ь — а)/К, существует и единствен периодический сплайн ¿т, м [У](х).
Доказательство. Рассмотрим систему (7). Умножим первую строку этой системы на матрицу и и сложим со второй строкой, полученную вторую строку умножим на матрицу и и сложим с третьей и т.д. Поменяем знаки. Тогда получим следующую систему:
/Е 0 . . 0 . .. —и ( X0 \ ( °0 \
0 Е. . 0 . .. —и2 X 2 А
0 0 . .Е. . —и1 X 2(г-1) = Д-1
\ 0 0 . . 0 . . Е — и V и 2^-1)) \Dl-j
(10)
1-1
где Д0 = В2А-1Рь-1, Д = V и) В2А-1 Р1-1-) — и1В2А-1 РЬ-1, I = 1,...,Ь — 1. По
)=0
условию теоремы
ёе^Е — иь) = 0 и X2(L-1) = (Е — иь)-1ВЬ-1, X21 = Д + и1+1 X2(L-1), и из (8) получаем X2г+1 = А-1(Р1 — AlX21) при I = 0,...,Ь — 1. Тем самым все коэффициенты периодического ¿-сплайна найдены.
3. Сходимость ¿-сплайна. Пусть периодическая функция /(х) € С6[а, Ь] и \/(хк) — Ук\ ^ ОН6+£, е ^ 0. Обозначим: <р{х) = ¡{х) - 5т>м(ж), Рг = <£>(&), <7г = ^'(б), П = \Ь? $"{£,{), г(х) = (<р(х), 1г<р'(х), ^Л2^ (ж)), XI = Приведем несколько вспомогательных утверждений, необходимых для доказатель-
ства оценок сходимости.
Лемма 1. Имеют место следующие соотношения:
г(х) = А(х)х1 + ш(х) + в(х)
х ^ С, I = 0,1,...,Ь — 1, где
, л (х — б)3 (х — б)4 (х — б)5
«11 (.Ж) — 1--—-¿378 --~л-¿638 --Г5-¿673,
Н3
Н4
Н5
— -1---^Я-¿478 --П-¿648 --^-¿674,
Н
Н3
Н4
Н5
, (х — Сг )2 (х — Сг )3 (х — Сг)4 (х — & )5
— -ТГГо---Гя-¿578 --~л-¿658 --Г^-¿675,
2Н2
Н3
Н4
Н5
1
а2г = ¡гаи, а:ц = |Н2аи (г = 1,2,3),
1
= ¿((ж - б)6 - Ь?(х - 6)^978 - Ь2{х - 6)^698 - Кх - 679),
,и)2(х) = киз^х), гУз(ж) = ^}г2'и)1(х).
Функции вг(х)/Н6 — 0 при Н — 0 для С ^ х ^ Сг + МН, г = 1, 2, 3. Доказательство. Воспользуемся формулой Тейлора для х ^ Сг:
Ф) =ш + + + (¿/'"(6) - 4) - б)3 + (¿/(4)(б) - 4) - б)4+
+ (^/(5)(6) - 4) (ж - б)5 + - б)6 + 0{(х- б)6) .
8 ВМУ, математика, механика, № 5
Рассмотрим разности К\ = ^/"'(^г) — а^, Я2 = ¿/^(б) — а\ и ^з = ^т/^Чб) — а\- С учетом выражений (2)—(4) получим
Т) 1 ГГ Л I 1 о- \ 1
= (6)-а3 = ^/ (й) - йе1;(А2)Лз
с1 5*7 Б С2 Б8 Б<9 С13 Б1°
(12)
Из формулы (3) следует
м
м
с1! = £[(Ут+к - /(т1 + к))к2+з] + £[(/(т1 + Л) - а° - а^Ьк - 4Ь2к2)к2+3], ; = 1,2,3
к=°
к=°
По предположению ут+к - /(т1 + к) = о(Ь6). Разложив /(т1 + к) в ряд Тейлора в точке б и подставив
это разложение в выражение для с,, получим
м
с, = £
к=°
/(б) - а° + (/ &) - а\)Ьк +
Лб) 2!
3!
+
-(М;)6 +о((М;)7) к2+^
4! 5! 6!
Подставляя выражения (10) для ро,Ц1,Г и используя обозначение (3), запишем
4 = Р^ 2+.? + + + К
4!
Б6+, +
++ 0(/1б)) ^ = 2)3
5!
6!
Подставим выражения для с, в (12). Получим
й'=- ь
Р^378 + &478 + ^¿578 + Л3 ¿678 + ^ + ^
/ (5)б )
3!
4!
5!
¿878 +
+^^^¿978 +о(/*7)
1
= + ^¿478 + П^578) ~ ^-^--Ь о(Л,4).
/ (6)т978
6!
Здесь мы воспользовались тем, что ¿678 = 1, ¿778 = ¿878 = 0. Аналогично
П ¿638 ¿648 ¿658 , 2 1 (бК , ги3\ п ¿673 ¿674 ¿675 , /(6'(б) , . ,,2\
2 = --- -^ГФ - -^-П - к 6, ¿698 + ), Дз = " - Д5"® " ~ 679 + о{к ).
Подставляя полученные выражения для Кг в (11), придем к требуемому равенству для ^>(х). Аналогично устанавливаются и два других равенства.
Лемма 2. Ясли р| < Сф6, ^\ < С^Ь6, \п\ < С3Ь6, то \ф(х)\ < С4Ь6, (х)| < С5Ь6, \ч>"(ж)| < С6Ь6 для ж € [б,б + МЬ], где константы С4,С5, С6 не зависят от Ь.
Лемма 3. Имеет место следующее рекуррентное соотношение:
гг+1 = и + щ + о(Ь6).
Здесь иг, = а^,(тЬ) — элементы матрицы и — имеют вид (9), од = щ(тЬ). Матрица и постоянная, не зависящая от свойств функции /(ж). Лемма 4. Если выполнено
\/(хк) - Ук\ < СЬ6+£, е > 0, (13)
и, кроме того, собственные значения матрицы и таковы, что |Аг| < 1, г = 1, 2, 3, то для периодического сплайна справедливы оценки
р\ < С1Ь6, \ < С2Ь6, \п| < С3Ь6, (14)
где константы С1,С2,С3 не зависят от I и Ь.
Доказательство. Из леммы 3 следует равенство
L—1
ZL = ^ икг^-1-к + иL¿0.
к=0
Так как для периодического сплайна ZL = ¿0, то ¿0 = (Е — и( ^ иkWL-l-k ).
к=0
Обозначим р = \\и||. Тогда
(Е — и^-1\\ < 1 + \\UL\\ + ... + \\UL\Г + ... <
1 — Р1
L—1
Еи к WZL—1—k к=0
<
L—1
ик
к=0
С7шаХ|/(6)|/.^С7таХ^1/(6)|"6 ' М 1 — Р
Отсюда следуют оценки (14) для Р0,Я0,Г0. Так как правая часть неравенства не зависит от индекса I, то неравенство (14) верно для любого I. Лемма доказана.
Теорема 2. Пусть периодическая функция /(х) € С6[а, Ь] и пусть выполнены предположения (13). Пусть, кроме того, собственные числа матрицы и по модулю меньше единицы. Тогда периодический сплайн ¿т м(х) с узлами на равномерной сетке имеет дефект 3 (т.е. ¿^м(х) € С2[а, Ь]) и для х € [а,Ь] справедливы следующие оценки:
ЯР - —
< СрН6-Р,
(15)
р = 0,1,2,3,4, 5; х = Сг при р = 3,4, 5; <^(р) (Сг) = ^(р) (Сг + О).
Доказательство. Из лемм 2, 3, 4 следуют оценки (15) при р = 0,1, 2 для х € [Сг, Сг+1], I = 0,1,... , Ь—1. По построению функция ¿т м(х) дважды непрерывно дифференцируема в узлах Сг, поэтому оценки (15) имеют место при р = 0,1, 2 для всех х € [а, Ь]. Применяя последовательно формулы численного дифференцирования, например
1
Ь?
[^(р)(х + Н) — 2^(Р) (х) + р(х — Н)] = ^(р+2)(х) + О(Н2),
х = Сг, ^(р+2)(Сг) = ^(р+2)(Сг + 0) при х = Сг, I ^ Ь — 1, и учитывая, что оценки (15) при р = 1, 2 уже доказаны, получим оценки при р = 3, 4, а затем и при р = 5.
Аналогичные утверждения справедливы и в непериодическом случае.
4. Устойчивость дважды непрерывно дифференцируемого ¿-сплайна. Собственные числа матрицы и определяются из уравнения
—Л3 + ^ их2 — (Д12 + Д13 + Л23 )Л + и = 0,
(16)
где и = чц + и22 + и33, Л) =
4гг 4г)
, ёе! и =
411 412 и 13 «21 422 423 . 431 «32 «33
Подставляя выражения для элементов матрицы и из (9) в уравнение (16), получим
9
Л3 — 3А2 + 3Л — 1 + агтг = 0,
г=1
где
а1 = (3^^578 Л2 + 3£578 — 6£578 Л),
а2 = ( — ^¿658 Л + 3£478 Л2 + 6£658 — 3£478 + 6£658 Л2),
а3 = (10*675 + ¿378 + 4*648 Л2 + 4*378 Л — 20*675 Л — 8*648 + 4*648 Л + ¿378Л2 + 10*675 Л2), а4 = (5*674 Л2 + 10*674 Л + 6*648*578 + 6*478 ¿658Л + 3*638 + ¿638Л2 — 15*674 — 6*478*658 — —6*578*648 Л + 8*638 Л),
1
«5 = (13Î673 Л + 6Î673 + 15Î674 ¿578 — 15Î478¿675 — 15Î674 ¿578Л + 15Î478 ¿675Л + ¿673 Л2 + +3Î37st658 Л — 3t638 ¿578 Л + 3t378t658 — 3t6381578),
а.6 = ( —10t6481675 — 10t674t658 Л — 1638t478 Л — 8t6731578 + 10t674 ¿658 — 1378t648 + 1378 ¿648Л—
— 7t673t578 Л + 1638t478 + 10t6481675Л + 8t3781675 + 7t378 ¿675Л), а.7 = ( — 2t673 ¿478 Л + 2t378t674 Л + 4t638 ¿675 Л + 3t6731478 — 3t3781674 — 6t673 ¿658 + 6t6381675 — 4t673t658 Л), «8 = ( — 3t638 ¿674 + 3t673 ¿648 — t67зt648 Л + t638t674 Л),
а9 = ( t378t674 ¿658 — 1673 ¿648't578 + ¿3781648'Î675 — ¿638 ¿478't675 + t638t674 ¿578 + ¿673~Î478 ¿658)-
Из уравнения (17) видно, что все три корня Л^^ — 1 при m — 0. Пусть / = Л — 1, ç = m/M, £ = 1/M. Для достаточно больших значений M получена асимптотика для собственных значений Л^:
/л3 + [27/2 — 27/4£ + 675/8£2 — 9747/16£3 + 162/35ç (1 — £ + 7£2 )](/2+
+ [432/5 — 432/5(1 — 7£ + 49£2 — 343£3 )£ + (—540/7 + 12960/91(1 — 7£ + 49£2 — 343£3 )£—
—2430/91(1 — 1/2£ + 1/4£2 )£)ç]ç2 / + [—4320/7 + 8640/7(1 — 7£ + 49£2 )£]( 2+
+[252 — 378£ + 3213£2 — 60669/2£3 + (—1080 + 2160(1 — 7£ + 49£2 )£)(4](3 = 0.
При M > 15 можно пользоваться следующим более грубым уравнением:
/3 + (27/2 — 27/4£)Z/2 + [432/5 — 432/5£ + (—1350/13 + 2025/13£)ç]ç2 / + (252 — 378£)(3 = 0. (18)
При Z = 1/15 корни уравнения
/3 + 27/2Z/2 + 432/5ç2 / + 252(3 = 0 (19)
равны —0, 41; —0, 25 — i ■ 0, 35; —0, 25 + i ■ 0, 35, что с точностью 0,1 отличается от корней (18). Уравнение (19) однородное. Обозначим v = //Z. Корнями уравнения
v3 + 27/2v2 + 432/5v + 252 = 0
являются vi = —6,11, V2 = —3, 69 — i ■ 5, 25, v-3 = —3, 69 + i ■ 5, 25. Отсюда получаем, что при достаточно малых m и больших M Л| = 1 — 6,11£ + O((2), |Л2,31 = 1 — 7, 38Z + O(Z2) по модулю меньше единицы. Тем самым доказана следующая теорема.
Теорема 3. Пусть ( = m/M < Z*. Тогда при достаточно малых m и больших M собственные числа матрицы устойчивости по .модулю меньше единицы.
Это условие устойчивости S-сплайнов аналогично условию устойчивости для кубического случая [1]. Для случая малых значений M (при 3 ^ M ^ 20) в результате расчета были получены значения собственных чисел матрицы U . Оказалось, что при m/M ^ < 1 все собственные числа матрицы U меньше единицы. Некоторые наиболее интересные значения m и M, при которых достигаются наименьшие значения max Л | и аппроксимация S-сплайнами устойчива, представлены в таблице.
Собственные числа матрицы U
M m Ai Л2 Аз max Ai m/M
4 2 -0,008 -0,231 -0,131г -0,231 +0,131г 0,265 0,25
5 3 -0,005 0,0549 - 0, 201г 0,0549 + 0,201г 0,207 0,60
6 2 0, 0266 -0,285 - 0,129г -0,285 + 0,129г 0,312 0,33
6 3 -0,008 -0,263 -0,0463г -0,263 + 0,0463г 0,266 0,50
7 2 0, 0732 -0,167 -0, 305г -0,167 + 0, 305г 0,347 0,285
7 4 -0,0069 -0,0737 - 0,214г —0,0737 + 0,214г 0,226 0,571
7 6 0,00218 0,116 -0,207г 0,116 + 0,207г 0,237 0,857
8 4 -0,0079 -0,265 - 0,031г -0,265 + 0,031г 0,266 0,50
8 5 -0, 00403 0,101 - 0,178г 0,101 +0,178г 0,204 0,625
8 7 0,00180 -0, 0466 - 0, 229г -0,0466 + 0,229г 0,233 0,875
9 5 -0, 00734 -0,124 -0,201г -0,124 + 0, 201г 0,235 0,555
9 8 0,00134 -0,205 — 0,118г -0,205 + 0,118г 0,236 0,888
10 5 -0,0078 -0, 263 - 0, 0407г -0, 263 + 0, 0407г 0,266 0,50
10 6 -0,0055 0,0182 - 0, 213г 0,0182 — 0,213г 0,213 0,60
11 7 -0,00322 0,141 - 0,147г 0,141 +0,147г 0,203 0,636
Все аналитические преобразования и вычисления на ЭВМ были произведены автором с помощью интегрированного математического пакета Мар1е [6].
СПИСОК ЛИТЕРАТУРЫ
1. Силаев Д.А, Якушина Г.И. Приближение ¿-сплайнами гладких функций // Тр. Семинара им. И.Г. Петровского. Вып. 10. М.: Изд-во МГУ, 1984. 197-206.
2. Силаев Д.А. Дважды непрерывно дифференцируемые ¿-сплайны // Вестн. Моск. ун-та. Матем. Механ. 2007. № 2. 12-17.
3. Силаев Д.А., Коротаев Д.О. Решение краевых задач с помощью ¿-сплайна // Математика. Компьютер. Образование: Сб. науч. трудов. Т. 2 / Под ред. Г.Ю. Ризниченко. М.; Ижевск: НИЦ "Регулярная и хаотическая динамика", 2006. 85-104.
4. Силаев Д.А, Амилющенко А.В., Лукьянов А.И., Коротаев Д.О. Полулокальные сглаживающие сплайны класса C1 // Тр. Семинара им. И.Г. Петровского. Вып. 26. М.: Изд-во МГУ, 2007. 347-367.
5. Silaev D.A., Amiliyushenko A.V., Luk'janov A.I., Korotaev D.O. Semilocal smoothing spline of class C1 // J. Math. Sci. 2007. 143, N 4. 3401-3414.
6. Матросов А. Maple 6. Решение задач высшей математики и механики. СПб.: БХВ-Петербург, 2001.
Поступила в редакцию 10.12.2007
УДК 519.6 -
БЛОЧНЫЕ ПРЕДОБУСЛОВЛИВАТЕЛИ КЛАССА ILU ДЛЯ ЗАДАЧ ФИЛЬТРАЦИИ МНОГОКОМПОНЕНТНОЙ СМЕСИ В ПОРИСТОЙ СРЕДЕ
К. Ю. Богачев1, Я. В. Жабицкий2
Рассмотрены предобусловливатели класса ILU (ILU(0), ILU(1), ILUT) для итерационных методов решения несимметричных линейных систем с разреженными матрицами, возникающими при аппроксимации систем дифференциальных уравнений в частных производных для задач фильтрации многокомпонентной смеси в пористой среде. Предложены блочный вариант хранения матриц и блочный вариант построения ILU-разложения, что позволяет значительно ускорить работу итерационных алгоритмов с ILU-предобусловливателем на рассматриваемом классе матриц. Результаты подтверждаются численными экспериментами с использованием различных матриц, полученных при дискретизации реальных задач на моделях нефтяных месторождений Западной Сибири.
Ключевые слова: блочное ILU-разложение, разреженные матрицы, итерационный алгоритм, предобусловливатель.
ILU class preconditioned (ILU(0), ILU(1), ILUT) employed for iterative algorithms for asymmetric linear systems with sparse matrices are considered. Test matrices used in this study originate from discretization of systems of partial differential equations describing a multi-component fluid flow in porous media. New algorithms for block storage of matrices and block based ILU-factorization are described. This new integrated approach was tested on a wide range of matrices resulted from actual hydrodynamic simulations of oil fields in Western Siberia and had demonstrated significant reduction of computational time.
Key words: block ILU factorization, sparse matrices, iterative algorithm, preconditioner.
Для решения систем с разреженной матрицей, возникающих при аппроксимации систем дифференциальных уравнений в частных производных, используются различные форматы хранения матриц, например Sparse или MSR [1]. Однако при таком подходе совершенно не учитывается, что в результате аппроксимации системы уравнений структуры ненулевых элементов каждого уравнения системы для одного блока сетки практически полностью совпадают. Именно поэтому возникла идея формализовать и
1 Богачев Кирилл Юрьевич — канд. физ.-мат. наук, доцент каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: bogachev@mech.math. msu. su.
2Жабицкий Яков Вячеславович — асп. каф. вычислительной математики мех.-мат. ф-та МГУ, e-mail: jjv@fromru.com.
