CC BY
35
Вычислительные технологии
Том 17, № 1, 2012
Сравнительный анализ реализаций модификации Рона в методах дробления параметров
Д. Ю. Людвин, С. П. Шарый Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: lyudvin@ngs.ru, shary@ict.nsc.ru
Работа посвящена сравнению различных реализаций методов дробления параметров для онтимальнохх) внешних) оценивания множеств решений интервальных .линейных систем уравнений. Главное внимание уделено анализу модификации на основе методики Рона как наиболее сложной, трудозатратной, но и наиболее эффективной. Представлены результаты вычислительных экспериментов с различными способами организации вычислений, их анализ и обсуждение. Даны практические рекомендации но оптимизации вычислительных схем методов дробления параметров.
Ключевые слова: интервал, интервальный анализ, интервальная система линейных алх'ебраичееких уравнений, метод дробления параметров, модификация Рона.
Введение и постановка задачи
Рассматривается задача внешнего интервального оценивания объединенного множества решений интервальных систем линейных алгебраических уравнений (ИСЛАУ) вида
anx i + aux2 + ... + аыхп = bi, tt2lXi + a22X2 + ... + a2nXn = b2, < : : , : : ^
anixi + an2X2 + ... + annXn = bn
с интервальными коэффициентами a^ и интервальными правыми частями bi; i,j = 1, 2,... ,n, или, кратко,
Ax = b, (2)
где A = (aij) — интервальная (пхп)-матрица, b = (bj) — интервальный n-вектор (брус). Системы (1)-(2) представляются как семейства обычных точечных (пеинтервальпых) линейных систем Ax = b той же структуры с матрицами A G A и векторами b G b.
Объединенным множеством решений интервальной линейной системы уравнений будем называть множество
^(A, b) = {x G Rn | (3 A G A)(3 b G b)(Ax = b)}, (3)
образованное всевозможными решениями точечных систем Ax = b с A G A, b G b. Для интервальных систем уравнений существуют также другие множества решений, соответствующие тем или иным практическим ситуациям. Поскольку в настоящей работе они не рассматриваются, то далее (3) будем называть множеством решений.
Далее интервальная матрица А предполагается неособенной, т. е. содержащей только неособенные (невырожденные) точечные матрицы. Тогда множество решений Б (А, 6) системы (1)-(2) ограничено. Нас будет интересовать задача оптимального внешнего интервального оценивания этого множества решений:
(4)
Вычисление дня множества решений внешних покоординатных оценок с заданной точностью является ХР-трудной задачей 11 — 31. В силу ее труднорешаемости соответствующие численные методы являются экспоненциально трудными и но своей структуре близки переборным алгоритмам дискретной оптимизации,
В работе рассматривается класс эффективных вычислительных алгоритмов дня нахождения решений задачи (4), называемый .методами дробления параметров (или РРБ-методами1) |4-6|. Их основная идея — представить "внешнюю задачу" как оптимизационную и применить дня ее решения интервальные методы глобальной оптимизации. Задача поиска нижних оценок множества решений ИСЛАУ эквивалентна задаче глобальной оптимизации величины
| Ах = Ь}
как функции параметров А € А и Ь € 6,
Задачу поиска нижних оценок множества решений ИСЛАУ можно переформулировать в виде
(5)
При поиске внешних оценок множества решений ИСЛАУ будем учитывать, что max{xv | x G Б (A, b)} = — min{xv | x G Б (A, — b)}.
Метод дробления параметров заключается в последовательном улучшении оценки min{xv | x G Б (A, b)} посредством дробления системы Ax = b на две ИСЛАУ-потомка, получающиеся из этой системы рассечением на концы одного интервального элемента Ab
согласно которой точные значения min{xv | x G Б (A, b)} и max{xv | x G Б (A, b)}, v = 1, 2,..., и, достигаются на решениях точечных систем уравнений Ax = b таких, что матрица A и вектор b образованы концами интервальных элементов из A и b. Итерационная процедура улучшения оценки строится в соответствии с хорошо известным методом "ветвей и границ".
Результат теоремы Бека-Никеля был значительно усилен И. Роиом |7|, уточнившим множество вершин матрицы A G A и вектора b G b, на которых достигаются минимальное и максимальное значения компонент точек из множества S(A, b) реше-
Ax = b
Цель данной работы — реализация различных версий методов дробления параметров, в частности, с использованием модификации Рона, их апробация на тестовых примерах и сравнительный анализ полученных результатов.
1PPS Partitioning Parameter Set.
Найти интервальный вектор и С Жга, имеющий наименьшую возможную ширину и содержащий множество решений Б (А, 6) интервальной системы ур авнений Ах = 6.
Найти min{xv | x G б (A, b)}, v =1, 2,..., и,
либо как можно более точные оценки снизу.
1. Метод дробления параметров для интервальных линейных систем уравнений
Существует ряд интервальных методов (например, интервальный метод Гаусса, интервальный метод Гаусса — Зейделя, метод Кравчика и т.д.), позволяющих получить внешнюю интервальную оценку множества S(A, b). Данные методы позволяют вычислить интервальный вектор, гарантированно содержащий множество решений ИСЛАУ, по не обеспечивают его оптимальность.
Обозначим через End какой-нибудь фиксированный метод внешнего оценивания
(A, b) G
IRn — интервальный вектор, получаемый с помощью этого метода, т. е.
End (A, b) D S(A, b).
Пусть Y(A, b) — нижний конец v-й компоненты (и = 1, 2,..., п) внешней интервальной оценки множества решений, получаемой методом End, т. е. Y(A, b) := (End(A, b))v. Основа PPS-методов — адаптивное дробление области параметров интервальной систе-
A
b
зультат.
Теорема Бека — Никеля [8, 9], Для, любого v G {1, 2,..., п} точные покоординатны,е оценки точек из объединенного мпооюеетва решений — экстремальные значат,я, min{ xv | x G S(A, b) } и max{ xv | x G S(A, b) } — достигаются, на, решениях край-
Ax = b A b
Ab
Потребуем от базового метода End удовлетворения условию: оценка Y(A, b) монотонна по включению относительно матрицы A и вектора b, т, е. для всех A', A'' G IRnxn и b', b'' G IRn при A' С A'' и b' С b'' верно неравенство
Y(A', b') > Y(A'', b''). (6)
В дальнейшем будем обозначать ИСЛАУ-потомки, получающиеся из интервальной системы Qx = r рассечением на концы одного интервального элемента либо в матрице Q, либо в векторе r, через Q'x = r' и Q''x = r'', Здесь Q' и Q'' — матрицы, полученные из Q заменой элемента q^ на q и q^ соответственно, а г' и г" — векторы, полученные из г заменой элемента па г^ и г» соответственно.
Решив две интервальных "системы-потомка" Q'x = r и Q''x = r, в силу свойства (6) можно получить более точную оценку снизу для искомого min{ xv | x G S(Q, r)} в виде
min{ Y(Q', r), Y(Q'', r) }.
Доказательство данного утверждения имеется в |4|, Аналогичный эффект дает и рас-
r ri
Zi и Vi.
Процедуру улучшения оценки для min{ xv | x G S(Q, r)} посредством дробления параметров можно повторить по отношению к системам-потомкам Q'x = r' и Q''x = r'', затем снова разбить потомков и улучшить оценку и т.д. Данную процедуру можно оформить в соответствии с методом "ветвей и границ".
2. Модификация Рона в методе дробления параметров
Известно, что объединенное множество решений интервальной линейной алгебраической системы имеет следующую элегантную характеризацию Оеттли — Прагера,
Теорема Оеттли — Прагера [10], Точка x £ Rn принадлежит множ еству Б (A, b) решений системы Ax = b тогда, и только тогда, когда,
| (mid A)x — mid b| < rad A ■ |x| + rad b. (7)
Суть теоремы Оеттли — Прагера состоит в том, что множество решений Б (A, b)
x
веиство, будем называть экстремальными решениями неравенства Оеттли — Прагера, Обозначим через E множество n-векторов с компонентами ±1, Очевидно, что множество E имеет мощно сть 2п, Для заданных век торов а, т £ E определим матрицы Ta , ACTT = {aj} и вектop Ь° = {Ь^} следующим образом:
Ta = diag {ai,... , о^},
AaT = mid A — Тст ■ rad A ■ TT, Ьа = mid b + Ta ■ rad b.
Из этого следует, что
I 0>г:, • еСЛИ GiTj = -1, I Ь, - еСЛИ <jj = 1,
aj 1 1 b = U i W
1 tty, если cfiTj = 1, I о,. если er,; = — 1.
Теорема Рона об экстремальных решениях [11], Пусть (и х и) -матрица A неособенна и b — интервальный и-вектор. Тогда для каждого а £ E уравнение
mid A ■ x — ■ rad A ■ |x| = ba (9)
имеет единственное решение xa, принадлежащее Б (A, b), и справедливо равенство
conv Б (A, b) = conv { xa | а £ E}, (10)
где conv — выпуклая оболочка.
Заметим, что вследствие (8) матрица AaT и вектор ba образованы наборами кон-Ab
Ax = b A
экстремальные решения могут достигаться лишь на множестве 4n матриц AaT и ассоциированных с ними векторов ba, т.е.
min{ xv | x £ Б (A, b)} = min ((ACTT )-1bCT )v , (11)
max{ xv | x £ Б (A, b)} = max ((ACTT )-1bCT )v (12)
для каждого v = 1, 2,..., и.
Описанный в раздело 1 метод дробления параметров можно модифицировать (см, |4|).
A
и вектора 6, на которых достигаются оптимальные внешние оценки множества решений ИСЛАУ,
Из (11) и (12) следует, что в процессе дробления исходной ИСЛАУ необходимо отслеживать концы задействованных интервальных элементов матрицы системы и вектора правых частей. Для этого с каждой интервальной линейной системой фж = г мы связываем
1) вспомогательную целочисленную (п х п)-матрицу Ш = |и>у}, элементы которой равны ±1 или 0, такую, что
(13)
2) вспомогательные целочисленные п-векторы в = {вг} и £ = {¿о} с компонентами ±1 0
Ыгз = вг£,- (14)
для всех г, ] = 1, 2,...,пи
если Яц аг0,
если
если Яго
если гг =
если гг = Ьг,
если гг = 6г.
(15)
Матрицу Ш и векторы в, £ будем называть контрольными. Значения определяются через матрицу Ш и вектор в посредством соотношения (14), Контрольные векторы в и £ являются "приближениями" к векторам о и т, входящим в равенства (11) и (12), а контрольная матрица Ш = в £т является "приближением" к матрице (о тт), В начале работы алгоритма элементы Ш, в и £ инициализируем нулями, далее в процессе дробления исходной ИСЛАУ они перевычисляютея.
Возникающие в процессе дробления системы-потомки организуем в список С, состоящий из записей вида
(ф, г, Т(ф, г),Ш,в,£, У, ж), (16)
где ф — интервальная (п х п)-матрица, ф С А; г — интервальный п-вектор, г С 6; Т(ф, г) — нижний конец г/-й компоненты внешней интервальной оценки множества решений Б(ф, г) Ш, в, £ — контрольные матрица и векторы; У — обратная интервальная матрица такая, что У Э | Q € ф}; ж — интерваль ный п-вектор такой, что ж Э Б(ф, г).
Образующие список С записи упорядочим по возрастанию значений оценки Т(ф, г). Первую запись списка, а также соответствующую ИСЛАУ фж = г и наименьшую в списке оценку Т(ф, г) будем называть ведущими на данном шаге.
На каждом шаге алгоритма с учетом значений контрольных матрицы Ш и вектора в выполняется процедура дробления интервального элемента Н ведущей ИСЛАУ фж = г. Если Н есть из матрицы ф (гк го вектора правой части г), то в случае и>к1 = 0 (вк = 0) порождаются две системы-потомка ф'ж = г' и = г". Если и>к1 = ±1 (вк = ±1), то оставляем в рабочем списке только одного потомка (какого именно, зависит от знака элемента ил и вк).
Заметим, что в процедуре дробления отбрасываются лишь те системы, которые либо не принадлежат множеству точечных систем х = | а, т € £}, либо не содержат такие системы. Поэтому отбрасывание систем-потомков не нарушает свойство ведущих оценок Т(ф, г) приближать искомый шт{ х^ | х € 5(А, 6)} снизу (предполагается, что оценка Т(ф, г) монотонна по включению относительно матрицы системы и вектора правых частей).
После дробления ведущей ПСЛАУ необходимо вычислить контрольные матрицы и векторы дня систем-потомков. Пусть в результате дробления порождаются две системы-потомка: ф'х = г' и ф''х = г". Обозначим соответствующие им новые контрольные матрицы через Ш 'и Ш'', а пары новых контрольных вектор ов — через з', ¿'и з'', ¿''. Новые контрольные матрицы и векторы вычислим следующим образом:
1) присваиваем Ш'' Ш' Ш, з'' з' з, ¿'' ¿'
2) а — если рассеченным элементом был из матрицы ф, то присваиваем и>к 1 и и>к — 1; б" — если рассеченным элементом был гк го вектора г, то присваиваем 4 <--1 и 4' ^ 1;
3) если хотя бы один из объектов семейства (Ш',з',£') изменился, перевычисляем остальные два, используя соотношение (14); эту процедуру повторяем, пока изменения не прекратятся; аналогично поступаем с (Ш'', з'', ¿'').
Если ведущая система порождает одного потомка ф'х = г', то контрольные матрица и векторы остаются неизменными, поэтому присваиваем Ш' Ш, з' в, ¿'
Отметим одно замечательное свойство контрольной матрицы Ш, вытекающее из соотношения (14), В каждой ее (2 х 2)-подматрице любой элемент равен произведению остальных трех элементов
= , (17)
иги'2 и»2Л и»1Л , (1®)
и»1Л и»2Л и»и2 • (1®)
Данные равенства могут быть полезными для уточнения контрольной матрицы Ш, Например, пусть мы намереваемся рассечь ведущую интервальную систему фх = г по элементу в то время как соответствующий элемент матрицы Ш равен нулю. При этом согласно обычному правилу дробления необходимо породить две системы-потомка. Но логично сделать попытку доопределить и^, поискав в Ш (2 х 2)-подматрицы, имеющие ненулевыми все элементы за исключением Если такая подматрица найдется, то присваиваем элементу значение произведения остальных трех элементов.
Заметим, что запись (ф, г, Т(ф, г), Ш, з, ¿, ^, ж), удовлетворяющая на некотором шаге условию
Х(д, г) >ш,
при дальнейшем выполнении алгоритма уже никогда не станет ведущей, поэтому либо вообще не вносится в список, либо удаляется при чистке списка, которая проводится после изменения (уменьшения) параметра ш.
На каждом шаге алгоритма проверяем монотонную зависимость решения системы от ее коэффициентов в соответствии с процедурой, описанной в |4-6|. На основе соотношений (13), (15) изменяем с 0 на ±1 значения элементов контрольных матрицы Ш и вектора з, соответствующие тем элементам в ф и г, на которых была выявлена монотонность.
В целях обеспечения наиболее быстрой сходимости алгоритма рекомендуется рассекать ведущие ИСЛАУ по элементам, па которых достигается максимум величин
дж^ (Q, r)
dq
v
wid q
iji
дж^ (Q, r)
dri
■ wid ri, i,j G {1, 2,...,n}, (20)
т. e, максимум произведения модуля интервальной оценки производной решения па ширину соответствующего интервала.
Достаточным условием оптимальности вычисленной оценки Y(Q, r), т.е. ее совпадения с min{ xv | x G Б (A, 6)}, является выполнение "условия точности": Y(Q,r) = (Q-1r)v для всех Q G Мгахга и r G Мга.
Наконец, опишем механизм контроля точности ведущих оценок, порождаемых ал-
Y(Q, r) Qx = r
также и величины Y(midQ,midr). Очевидно, что
Y(mid Q, mid r) > Y(Q, r),
и значения Y(mid Q, mid ^приближают min{ xv | x G Б (A, 6)} сверху: если для каждого шага алгоритма определять
ш = minY(mid Q, mid r) (21)
Qx = r
L до этого шага, то
min{ xv | x G Б (A, 6)} < ш. Qx = r
Y(Q, r) < min{ xv | x G Б (A, 6)},
и потому одним из критериев остановки алгоритма может быть достижение требуемой малости величины (ш — Y(Q, r)).
Итоговая схема алгоритма метода дробления параметров с модификацией Рона приведена ниже (см. с. 76).
3. Реализация алгоритмов и численные эксперименты
Представленный в предыдущем раздело алгоритм был реализован в программной системе Matlab, дополненной пакетом интервального расширения IXTLAB 5.5. Помимо основной схемы алгоритма внешнего оценивания множества решений ИСЛАУ методом дробления параметров с модификацией Рона были реализованы его модификации, различающиеся базовым алгоритмом End, а также способом организации и обработки списка. Базовый алгоритм End реализован на основе
— метода Кравчика |12|,
- модифицированного метода Кравчика с использованием "эпеилои-разд\'тия" |13, 14|,
— интервального метода Гаусса |4, 15|,
— интервального метода Гаусса — Зейделя |4|,
— процедуры Хансена —Блика —Рона (изложена в |11, 16|),
Алгоритм оптимального внешнего оценивания множества решений ИСЛАУ методом дробления параметров с модификацией Рона
DO WHILE ((ведущая оценка Y(Q, r) неточна) и (w — Y(Q, r) > е))
дх (^^ r)
вычисляем интервальные расширения производных ——--— и
dqij
дх^(^Q r)
——■— но элементам q.r. и ri с ненулевой шириной; dri j
сжимаем в Q и r элементы, на которых выявлена монотонная
зависимость xv от qj и ri ; соответствующие элементы в контрольных матрице W и векторе s
изменяем с 0 на ±1 на основе (13), (15); ищем в ведущей ИСЛАУ Qx = r интервальный элемент h,
которому соответствует наибольшее из произведений (20);
W
порождаем одну или две ИСЛАУ-потомка Q'x = r' и Q''x = r''; если порождены две системы-потомка, неревычисляем матрицы W', W'' и векторы s', s'', t', t'' согласно (14), иначе присваиваем W' W, s' ^ s t' ^ t; вычисляем вектор ж' ^ Encl(Q', r') и, возможно, ж'' ^ Encl(Q", r''); присваиваем оценку v' ^ Y(Q', r') и, возможно, v'' Y(Q'', r''); вычисляем оценку обратной интервальной матрицы Y' 5 (Q')-1
и, возможно, Y'' 5 (Q'')-1; вычисляем оценку Y(mid Q', mid r') и, возможно, Y(mid Q'', mid r'') и присваиваем ^ ^ min{ Y(mid Q', mid r'), Y(mid Q'', mid r'')}; удаляем из списка L бывшую ведущую запись (Q, r, v, W, s, t, Y, ж); если v' < w, т0 помещаем запись (Q', r', v', W', s', t', Y', ж') в список L
v
если ведущая ИСЛАУ породила две системы-потомка и v'' < w,
то помещаем запись (Q'', r'', v'', W'', s'', t'', Y'', ж'') в список L
v
если w > то полагаем w ^ ^ и чистим спи сок L: удаляем из него все такие записи (Q, r, v, W, s, t, Y, ж), что v > w;
END DO
В качестве алгоритма End была использована также процедура verif ylss из пакета INTLAB [17], Процедура x = verif ylss (A, b) предназначена для внешнего оценивания множества решений интервальной системы Ax = b и состоит из двух этапов, первый из которых реализует алгоритм, основанный на операторе Кравчика, второй — процедуру
x
в течение семи итераций не находится, то выполняется второй этан, т. е, решение будет получено с помощью процедуры Хансена —Блика —Рона.
Следует отметить, что во всех базовых алгоритмах применялось предобуславлива-ние ИСЛАУ обратной средней матрицей. Соответствующие базовые алгоритмы использовались также дня внешнего оценивания "обратной интервальной матрицы", которая применяется в тесте на монотонность и при выборе рассекаемого интервального элемента.
Реализованы следующие способы организации и обработки списка:
- список С упорядочен по возрастанию оценки Т(ф, г) таким образом, что первая его запись является одновременно и ведущей. Чистка бесперспективных записей производится только при изменении параметра ш;
- список С организован в виде кучи [18]. Ведущая запись, для которой оценка Т(ф, г) минимальна, определяется при просмотре всего списка. Чистка бесперспективных записей осуществляется только при изменении параметра ш;
- в списке С выделяется упорядоченный подсписок С7 = {(ф, г, Т(ф, г),Ш,
У, ж) | Т(ф, г) < 7} записей, которые будем называть активными. Здесь 7 — некоторая "пороговая константа" такая, что Т(ф, г) < 7 < ш, Дополнение С\С7 организовано в виде кучи. Если в процессе работы алгоритма подмножество С7 станет пустым, то по-
роговая константа 7 перевычисляется. Например, можно взять 7 =
1
1 + А
(А-ш+Т(ф, г)),
где 0 < А < 1, и из С снова выделить упорядоченный подсписок С7, Производится чистка пе всего списка, а только подсписка активных записей. Данный способ обработки списка предложен П.С. Пайковым |19|;
- в списке С выделяется упорядоченный подсписок Сг активных записей, имеющий некоторую фиксированную максимальную длину /; дополнение С = С\Сг организовано в виде кучи. Ведущей записью считается первая запись Сг,
Если на некотором шаге алгоритма длина Сг станет равной /, т. е. подсписок актив-пых записей полностью заполнится, то новая запись в зависимости от содержащейся в ней оценки Т(ф, г) может быть внесена либо в Сг, либо в подсписок С, организованный в виде кучи. Обозначим через Ттах максимальную оценку Т(ф, г), содержащуюся в записях подсписка Сг, В случае, когда длина Сг равна /, новая запись с оценкой Т(ф, г) < Ттах вносится в подсписок активных записей. При этом происходит его переполнение и лишняя запись помещается в кучу С. Записи-потомки с оценками Т(ф, г) > Ттах при полностью заполненном подсписке Сг вносятся в С.
Если подсписок Сг активных записей становится пустым, то из С отбираются не более / записей с наименьшими оценками Т(ф, г) и помещаются упорядоченно в Сг,
Чистка подсписка С производится при изменении п араметра ш, чистка подсп иска Сг выполняется только при С = 0,
Представленные алгоритмы были апробированы па следующих тестовых интервальных линейных системах уравнений.
Пример 1. Интервальная линейная система Ноймайера |20| имеет вид
( 0 [0,2] [0,2] 0
\[0, 2] [0, 2]
[0, 2] \
[0, 2]
0 )
х
/[—1,1]\
[—1,1]
1,1]/
(22)
где 9 — неотрицательный вещественный параметр. Отметим, что матрица системы четного порядка п пеособеппа при 0 > п, нечетного порядка п — при 0 > у/п2 — 1. При
приближении 9 к границам неособенности размеры объединенного множества решений ИСЛАУ (22) неограниченно возрастают. Варьируя 9, можно получить набор тестов для проверки рассмотренных алгоритмов оптимального решения "внешней задачи" ИСЛАУ, Пример 2. Интервальная линейная система уравнений, впервые рассмотренная С,И, Шарым в |4|, имеет вид
( [п - 1, N [а - 1,1 - в] [а - 1,1 - в] [п - 1, N]
[а - 1,1 - в] \ [а - 1,1 - в]
х
\ [а - 1,1 - в] [а - 1,1 - в] ... [п - 1, N ]
( [1 - п,п - 1] \ [1 - п, п - 1]
V [1 - п,п - 1]
т
где п — размерность системы (п > 2), 0 < а < в < 1 N — вещественное число, не меньшее п - 1, Варьируя значения а., в, п и N, получим широкий набор ИСЛАУ
в
пуню, матрица системы (23) становится все более близкой к особенной, а множество
ав
модифицировать форму множества решений.
Можно показать [4], что оптимальными покомпонентными оценками множества Б решений системы (23) являются
тт{ хг | х € Б } = -1/а,
тах{ хг | х € б } = 1/а, г = 1, 2,
п,
причем эти оценки не зависят от конкретного значения N.
Пример 3. Интервальная линейная система Тофта |21| имеет вид
Ах = Ь,
(24)
где матрица системы
( [1 - г, 1 + г] 0
0 [1 - г, 1 + г]
00
\ [1 - г, 1 + г] [2 - г, 2 + г]
0 0
[1 - г, 1 + г] [2 - г, 2 + г]
[1 - г, 1 + г] [п - 1 - г, п - 1 + г] [п - 1 - г, п - 1 + г] [п - г, п + г] У
и вектор правой части
( [1 - Я, 1 + Я] \ [1 - Я, 1 + Я]
V [1 - я, 1 + Я] у
Здесь г Я — положительные вещественные числа. Матрица данной системы является интервализацией известной тестовой матрицы вычислительной линейной алгебры из справочника |22|,
4. Сравнительный анализ результатов 4.1. Влияние базового алгоритма
Мы провели ряд тестовых расчетов с интервальной системой Ноймайера (см, пример 1), меняя размерность системы и, значение ее параметра в и базовый алгоритм End внешней оценки множества решений интервальной системы, В табл. 1 представлены результаты численных экспериментов — время нахождения оптимальной нижней оценки первой компоненты множества решений ИСЛАУ (см. пример 1). Базовые алгоритмы обозначены следующим образом: К — метод Кравчика, МК — модифицированный метод Кравчика, G — интервальный метод Гаусса, GS — интервальный метод Гаусса — Зейделя, HBR — процедура Хансена —Блика —Рона, V — процедура verifylss из пакета IXTLAB.
Проанализировав полученные результаты, можно сделать следующие выводы.
Метод Кравчика, являясь итерационным, в качестве базового алгоритма оказался наименее эффективным. Несколько быстрее работает программа, если базовый алгоритм основан па модификации метода Кравчика и па интервальных методах Гаусса
Т а б .л и ц а 1. Сравнение базовых алгоритмов (ем. пример 1)
Параметр в
Время работы программы на основе базовых алгоритмов, с
матрицы системы К МК G GS HBR V
п = 5
10 3.314 2.847 2.814 2.818 1.046 2.009
14 0.386 0.326 0.700 0.714 0.284 0.338
18 0.269 0.259 0.569 0.614 0.221 0.270
22 0.210 0.229 0.500 0.571 0.195 0.238
26 0.212 0.199 0.499 0.539 0.189 0.207
30 0.214 0.197 0.499 0.547 0.188 0.204
п = 8
16 390.255 103.098 303.903 289.908 65.099 89.632
20 26.002 7.016 17.196 17.372 5.432 6.144
24 3.811 2.604 6.305 6.544 2.004 1.881
28 1.977 1.744 4.730 4.817 1.500 1.408
32 1.384 1.357 3.479 3.578 1.104 1.038
36 1.088 1.301 2.736 2.791 0.838 0.789
40 0.894 0.962 2.265 2.283 0.688 0.649
44 0.844 0.813 2.107 2.091 0.641 0.604
48 0.743 0.763 1.948 1.971 0.589 0.557
п = 10
22 3282.707 610.170 1142.017 1028.602 166.874 520.371
26 100.611 30.148 84.395 82.148 25.425 26.211
30 20.640 12.953 34.359 34.145 10.293 9.4704
34 11.968 6.359 18.369 18.371 5.392 4.9562
38 6.561 4.994 14.252 14.230 4.185 3.852
42 5.435 3.719 11.098 11.075 3.256 3.003
46 3.865 3.182 8.694 8.661 2.555 2.423
50 2.836 2.595 8.459 8.408 2.474 2.290
54 2.681 2.516 7.017 6.969 2.058 1.898
58 2.286 2.311 5.813 5.750 1.704 1.573
и Гаусса — Зейделя, Применение этих алгоритмов как базовых, скорее всего, нежелательно.
Наилучшую скорость работы показала программа, использующая процедуру Хансена— Блика — Рона. Эффективность применения в качестве базового алгоритма процедуры verifylss из пакета IXTLAB оказалась несколько хуже за счет того, что па первом этане ее работы, основанном па методе Кравчика, решение не достигается, после чего используется процедура Хансена — Блика — Рона,
Таким образом, при выборе базового алгоритма следует предпочесть процедуру Хансена — Блика — Рона.
При анализе примера 1, изменяя значение параметра в, было исследовано влияние свойств интервальной матрицы системы — ее близость к границам неособенности -на работу программы. В табл. 1 приведены результаты численных экспериментов с
в
те из них, при которых время работы программы не слишком велико).
Дня характеристики свойств интервальной матрицы системы использованы вони-чины:
1) р = р (| (mid A)-1| • rad A) — спектральный радиус матрицы |(mid A)-1| • rad A;
2) Да = am¡n(mid A) — amax(rad A) — разность между наименьшим сингулярным
mid A rad A
Напомним, что сингулярными числами матрицы A е Rnxn называются неотрицательные квадратные корни из собственных чисел матрицы AAT,
Выбор указанных выше величии мотивируется следующими признаками неособенности интервальных матриц.
Признак Риса — Бека. Пусть интервальная (n х n)-матрица A такова, что mid A неособенная и р (|(mid A)-1| • rad A) < 1. Тогдa A неособенная.
(n х n) A mid A
и max (rad A • |(mid A)-1|).. > 1. Тогда, A особенная.
1<j<n jj
Признак Румпа, Пусть для, интервальной (пхп)-матрицы A имеет место неравенство amax(rad A) < amin(mid A). Тогд a A неособенная.
Признак Рона — Рекса. Если, для, интервальной (n х п)-жатрицы A имеет место amin(rad A) > amax(mid A), то она особенная.
В табл. 2 приведены значения рассматриваемых характеристик р и Да матрицы Ноймайера при различных значениях ее диагонального параметра в и размерности п.
Отметим, что при близких к нулю значениях р и достаточно больших Да время работы программы мало. Однако оно экспоненциально возрастает при приближении матрицы системы к границам неособенности, т.е. при р — 1, Да — 0.
р
и Да интервальной матрицы A размерности п = 5 для различных базовых алгоритмов.
В табл. 3 представлены результаты численных экспериментов с интервальной си-
n
ее параметров а в N и базовых методов внешней оценки множества решений интер-
р Да
системы и время нахождения оптимальной нижней оценки первой компоненты множества решений ИСЛАУ на основе разных базовых алгоритмов. На рис. 2 изображены диаграммы эффективности работы модификаций программы при разных значениях параметров и размерности интервальной системы Шарого,
Таблица 2. Характеристики р и А а матрицы Ноймайера
Параметр в р(|(гшс1А) 1 -гаёА) <7тт(1шё А) - атах(гаё А)
п = 5
10 0.5397 5
14 0.3590 9
18 0.2674 13
22 0.2125 17
26 0.1760 21
30 0.1501 25
п = 8
16 0.5884 8
20 0.4503 12
24 0.3633 16
28 0.3037 20
32 0.2605 24
36 0.2279 28
40 0.2024 32
44 0.1819 36
48 0.1652 40
п = 10
24 0.5392 14
28 0.4423 18
32 0.3740 22
36 0.3235 26
40 0.2846 30
44 0.2539 34
48 0.2291 38
52 0.2086 42
56 0.1914 46
60 а 0.1767 50 5
0.3 р 0.4
Рис. 1. Зависимость времени работы программы £ от характеристик р (а) и Аа (б) матрицы Ноймайера
Анализируя полученные результаты, можно сделать ьыьоды, аналогичные изложенным выше. Наименее эффективной оказалась программа, использующая в качестве базового метод Кравчика и его модификацию (эти данные в табл. 3 не приводятся). Наилучшую скорость работы имеем в случае, когда базовым методом является процедура Хансена — Блика — Рона. Скорость работы программы па основе процедуры уег1£у1ББ из пакета ТХТЬАВ значительно надает при возрастании размерности интервальной си-
Т а б .л и ц а 3. Сравнение базовых алх'оритмов (пример 2)
Параметр N Характеристика Время работы программы
интервальной матрицы системы иа основе базовых алх'оритмов, с
матрицы системы Р До- НВИ V
а = 0.4, в = 0.6, п = 10
15 0.6757 3.6 11.225 9.054 2.356 3.578
20 0.7353 3.6 10.823 9.322 2.341 3.853
25 0.7764 3.6 10.811 9.211 2.322 3.855
а = 0.4, в = 0.6, п = 20
25 0.6219 7.6 105.509 93.826 21.935 39.846
30 0.6637 7.6 105.558 97.270 21.959 25.755
35 0.6972 7.6 105.554 95.336 21.880 77.288
а = 0.4, в = 0.6, п = 30
35 0.6014 11.6 386.129 407.671 91.435 456.992
40 0.6329 11.6 386.522 426.977 91.432 358.437
45 0.6598 11.6 385.823 430.413 92.363 676.331
а = 0.6, в = 0.8, п = 10
15 0.5150 5.4 3.659 3.582 1.081 2.313
20 0.6029 5.4 7.915 6.450 1.622 1.876
25 0.6646 5.4 7.904 6.377 1.626 1.919
а = 0.6, в = 0.8, п = 20
25 0.4328 11.4 23.631 21.467 5.309 17.544
30 0.4956 11.4 23.892 21.540 5.310 47.364
35 0.5458 11.4 27.784 24.573 6.299 97.143
а = 0.6, в = 0.8, п = 30
35 0.4021 17.4 74.710 68.064 16.188 216.368
40 0.4494 17.4 74.648 67.998 16.192 441.926
45 0.4897 17.4 74.671 67.643 16.188 1616.238
а = 0.3, в = 0.7, п = 10
15 0.7353 2.7 13.542 12.483 3.373 5.410
20 0.7874 2.7 13.435 12.864 3.259 4.170
25 0.8224 2.7 13.457 12.525 3.360 4.186
а = 0.3, в = 0.7, п = 20
25 0.6868 5.7 157.602 151.104 34.394 34.016
30 0.7246 5.7 159.019 144.843 34.416 34.065
35 0.7543 5.7 161.887 141.492 34.346 37.167
а = 0.3, в = 0.7, п = 30
35 0.6679 8.7 712.694 656.571 149.929 145.106
40 0.6969 8.7 704.724 670.453 150.615 146.732
45 0.7212 8.7 712.297 647.965 151.154 147.519
t, С
t, с 120
Л/=25 Л/=30 Л/=35
и = 20, 0.6, Р= 0.8
Л/=25 Л/=30 Л/=35
и = 20, ог= 0.4, Р= 0.6
Л/=35 Л/=40 Л/=45
и = 30, а= 0.6, Р= 0.8
Л/=35 Л/=40 Л/=45
п = 30, ог= 0.4,/?= 0.6
Рис. 2. Эффективность а;п'оритма при разных значениях параметров и размерности интервальной системы Шарохх)
егемы за счет увеличения времени выполнения итераций метода Кравчика на первом этапе. Модификации на основе базовых методов Гаусса и Гаусса — Зейделя показывают достаточно близкие по эффективности результаты.
Из табл. 3 также видно, что при фиксированных значениях и, а и в увеличение параметра N не приводит к изменению характеристики Да интервальной матрицы системы (при этом р варьирует в небольших пределах). Поэтому скорость работы программы, основанной па базовых методах Гаусса, Гаусса — Зейделя и процедуре Хансена — Блика — Рона, меняется незначительно. Однако модификация программы, испо.иь-
N
рис. 2). Это объясняется тем, что интервальная оценка решения, как правило, вычисляется па первом этапе работы процедуры, т. е. с помощью модифицированного метода
N
Далее была проведена апробация разработанных алгоритмов и их модификаций па тестовой интервальной системе, описанной в примере 3 (табл. 4). Сравнивая результаты работы алгоритмов, основанных па разных базовых методах, можно сделать выводы, аналогичные полученным рапсе.
Таким образом, па ряде тестовых примеров апробированы шесть базовых методов дня использования их в разработанном алгоритме дробления параметров с модифика-
Т а б .л и ц а 4. Сравнение базовых алх'оритмов (пример 3)
Параметры Характеристика Время работы программы
матрицы системы матрицы системы иа основе базовых алх'оритмов, с
г = К Р в НВИ V
п = 5
0.1 0.2281 0.480 0.311 0.221 0.378
0.2 0.4562 0.422 0.375 0.116 0.199
0.3 0.6844 1.161 1.151 0.307 0.688
п = 10
0.1 0.2028 0.813 0.810 0.200 0.216
0.2 0.4056 2.185 1.921 0.475 0.435
0.3 0.6083 4.083 5.345 0.820 1.399
п = 15
0.1 0.1998 1.771 1.769 0.419 0.410
0.2 0.3997 5.861 5.101 1.154 0.838
0.3 0.5995 13.834 16.335 2.594 4.690
п = 20
0.1 0.1991 3.107 3.094 0.696 0.678
0.2 0.3983 9.126 7.789 1.571 1.728
0.3 0.5974 41.162 43.119 7.576 11.403
п = 25
0.1 0.199 4.863 4.805 1.062 1.592
0.2 0.3979 14.129 12.084 1.694 2.016
0.3 0.5969 90.467 104.553 16.273 27.129
п = 30
0.1 0.1989 10.969 10.546 1.557 2.238
0.2 0.3979 22.807 19.828 3.450 4.110
0.3 0.5968 216.787 231.572 39.057 60.433
дней Рона, Среди них наиболее эффективной оказалась процедура Хансена —Блика — Рона, Время работы алгоритма нахождения оптимальной внешней оценки множества решений ИСЛАУ зависит как от ее размерности, так и от близости к границам неособе шюети,
4.2. Эффективность процедуры работы со списком
В раздело 3 были описаны следующие четыре способа организации и обработки списка записей, в которых хранится информация о системах-потомках, полученных в результате дробления параметров:
1) список С организован в виде кучи;
2) записи списка С упорядочены по возрасталию оценки Т^, г);
3) в списке С выделяется упорядоченный подсписок активных записей, имеющий некоторую фиксированную максимальную длину, а остальные записи организованы в виде кучи;
4) способ, предложенный Пайковым,
Соответствующие модификации алгоритма оитималыюго внешнего оценивания множества решений ИСЛАУ, реализующие данные способы работы со списком, аиробиро-
Т а б .л и ц а 5. Эффективность алгоритмов организации и обработки списка
Параметры и размерность Время работы программы, с
интервальной системы способ 1 способ 2 способ 3 способ 4
Интервальная система Нойлшйе.ра,
п = 6 в = 12 4.347 3.013 2.082 2.943
в = 16 1.000 0.677 0.667 0.670
в = 20 0.596 0.404 0.399 0.400
п = 8 9 = 16 119.767 73.602 65.039 72.844
в = 20 10.241 5.452 5.432 5.427
в = 24 3.859 2.036 2.004 2.031
п = 10 9 = 20 5815.058 2323.042 1407.043 2235.154
в = 22 336.505 185.925 166.874 184.065
в = 24 67.605 41.349 40.886 41.154
п = 12 в = 28 853.403 461.624 420.766 460.146
в = 30 301.189 178.815 175.215 176.626
в = 35 84.595 48.643 48.003 48.689
Интервальная система Шарого
п = 6 N = 20, п = 10 5.739 3.428 3.259 3.275
N = 30, п = 20 67.014 35.038 34.416 35.523
N = 40, п = 30 297.970 151.301 150.615 151.303
Интервальная система Тофта
п = 8 п = 10 1.340 0.856 0.820 0.857
п = 20 13.358 7.772 7.576 7.784
п = 30 71.625 39.921 39.057 39.930
вапы на тестовых системах (см, примеры 1-3). Результаты численных экспериментов представлены в табл. 5. В качестве базового алгоритма использовалась процедура Хансена — Блика — Рона.
Наименее эффективным оказался способ организации списка в виде кучи (способ 1). Очевидно, па каждом шаге алгоритма достаточно много времени требуется па просмотр всего списка с цслыо поиска ведущей записи. Скорость обработки списка значительно повышается, если список упорядочивается по возрастанию оценки Т^, г) (способ 2). Однако процесс сортировки списка, добавления и удаления записей из пего требует определенных временных затрат. Поэтому имеет смысл поддерживать упорядоченность, а также производить чистку не всего списка С, а некоторой его части -подсписка активных записей — способы 3 и 4. Если последние сравнивать по временным затратам, то следует отдать предпочтение способу 3, при котором фиксируется максимальная длина подсписка активных записей. В способе 4, предложенном Папко-вым, длина подсписка активных записей определяется пороговой константой 7 и на некоторых шагах алгоритма может быть достаточно большой, что приводит к снижению его быстродействия, особенно в случае, когда матрица системы близка к границам пеособеппости.
При организации списка по способу 4 па каждом шаге алгоритма ведущая запись содержит наименьшую оценку Т^, г). Данное свойство может нарушаться при организации списка по способу 3, что приводит к увеличению количества итераций, необходимых дня получения оптимального решения ПСЛАУ,
С
900
800
700
■ 1 600
2 500
■ 3 400
■ 4 300
200
100
0
Рис. 3. Эффективность алгоритмов организации и обработки списка
По результатам тестовых экспериментов с системой Ноймайера при разных размерностях (п = 6 и п = 12) и значениях пара метра 0 построены диаграммы (рис. 3). Напомним, что при возрастании 0 увеличивается характеристика матрицы системы р и уменьшается характеристика Да (см. табл. 2), т.е. матрица ИСЛАУ приближается к границам неособенности. При этом, как видно из рис. 3, разница в быстродействии различных способов обработки списка становится более заметной.
4.3. Сравнение процедур оптимального внешнего оценивания множества решений ИСЛАУ
Приведем результаты сравнительного анализа работы двух алгоритмов оптимального внешнего оценивания множества решений интервальной линейной системы уравнений Ах = Ь:
- процедуры уегз.1гЬегуа11ш11 из пакета УЕИБОЕТ |23|, основанной па методе Рона (см. раздел 2) и написанной им самим;
- процедуры ПпррБг, реализующей метод дробления параметров с модификацией Рона, где в качестве базового метода используется процедура Хансена — Блика — Рона, а список записей организован и обрабатывается по способу 3 (см. раздел 4.2).
В таб.и. 6 представлены результаты тестовых расчетов с интервальными системами (см. примеры 1-3) при различных размерностях и значениях параметров матрицы ИСЛАУ. Следует заметить, что задача вычисления оптимальных внешних оценок множества решений интервальной системы является ХР-трудной. Время работы рассмот-
п
мы. Из таблицы видно, что скорость работы процедуры уегз.1гЬегуа11ш11 не зависит
А
различных значениях ее параметров. Быстродействие алгоритма Ппррвг снижается, А
п
Как было показано в разделе 2, искомые оценки тш{ х„ | х € 5(А, Ь)} и тах{ х„ | х € 5(А, Ь)}, V = 1, 2,..., п, достигаются на множестве не более чем 2п экстремальных решений, т. е. верхняя оценка сложности метода Рона равна 2п, Сложность методов
2п
Т а б .л и ц а 6. Эффективность алх'оритмов verintervalhull и linppsr
Параметры и размерность интервальной системы Время работы алх'оритма, с
linppsr verintervalhull
Интервальная система Ноим.аиера,
п = 5 в = 10 8.012 3.017
в = 20 1.642 3.019
в = 30 1.641 3.025
п = 10 в = 25 617.238 144.113
в = 30 197.024 144.226
9 = 45 51.422 144.815
в = 60 31.065 145.012
п = 12 в = 30 2943.612 1084.556
в = 50 246.804 1084.476
в = 70 99.308 1085.555
в = 90 57.049 1084.001
Интервальная система Шарого
а = 0.4, ¡5 = 0.6 п = 10, N = 15 42.528 144.717
и = 20, N = 25 882.344 Более 8 ч
п = 30, N = 35 5483.106 Более 8 ч
а = 0.6, ¡5 = 0.8 п = 10, N = 15 17.477 142.084
и = 20, N = 25 209.250 Более 8 ч
и = 30, N = 35 969.280 Более 8 ч
Интервальная система Тофта
г = 0.2 п = 10 5.750 9.115
п = 20 56.524 1321.249
п = 30 161.393 Более 8 ч
Однако менее трудоемким алгоритм дробления параметров становится за счет того, что па каждом его шаге существенно используется информация, полученная в ходе выполнения предыдущих шагов, т.е. алгоритм в данном случае является адаптивным.
Важная особенность метода дробления параметров состоит также в том, что оп порождает последовательность приближенных оценок для min{ xv | x £ £ (A, b)} и max{ xv | x £ £ (A, b)}, v = 1, 2,... ,n, снизу и сверху соответственно. Поэтому в случае трудпорешаемой задачи, когда размерность системы достаточно велика или матрица системы близка к границам пеособеппости, процесс работы алгоритма linppsr может быть прерван до своего полного завершения, причем полученные оценки будучи по оптимальными могут служить более или менее точными приближениями искомых min{ xv | x £ £ (A, b)} и max{ xv | x £ £ (A, b)}, v = 1, 2,..., n.
Таким образом, среди рассмотренных реализаций методов дробления параметров с модификацией Pona, различающихся базовым алгоритмом (метод Кравчика, модифицированный метод Кравчика с использованием "эпеилоп-раздутия", интервальный метод Гаусса, интервальный метод Гаусса — Зейделя, процедура Хансена —Блика —Рона), наиболее предпочтителен по времени работы алгоритм, использующий процедуру Хап-сепа — Блика — Pona с предобуславливапием обратной средней матрицей. Дня достижения наибольшей эффективности алгоритма имеет смысл при работе со списком систем-потомков, полученных в результате дробления параметров, выделять подсписок ак-
тивпых записей, в котором поддерживается упорядоченность по значениям оценки и производится чистка от "бесперспективных" потомков.
Анализ результатов апробации алгоритма linppsr, реализующего метод дробления параметров с модификацией Рона, показал, что время работы программы существенно зависит от свойств интервальной матрицы системы, в частности, от близости этой матрицы к границе множества неособенных (невырожденных) матриц. Если матрица системы "почти особенная", а также если размерность системы велика, быстродействие алгоритма снижается, В этом случае дня получения решения за практически приемлемое время можно воспользоваться свойством последовательной гарантии методов дробления параметров и прервать работу программы дня вычисления некоторого приближения искомой внешней оценки. Это свойство алгоритма linppsr выгодно отличает его от других подходов к оптимальному внешнему оцениванию множеств решений интервальных линейных систем, например, рассмотренной выше процедуры verintervalhnll, реализующей метод Рона.
Список литературы
[1] Kreixovicii V., Lakeyev A.v., Roiix J., Kahl P. Computational Complexity and Feasibility of Data Processing and Interval Computations. Dordrecht: Kluwer, 1997.
[2] KREixovicii V., Lakeyev A.V., Noskov S.I. Optimal solution of interval linear systems is intractable (NP-hard) /7 Interval Comput. 1993. No. 1. P. 6 14.
[3] Ronx J., KREixovicii V. Computing exact componentwise bounds on solutions of linear systems with interval data is NP-hard /7 SIAM J. on Matrix Analysis and Appl. 1995. Vol. 16. P. 415 420.
[4| Шарый С.П. Конечномерный интервальный анализ. Электронная книга http://www-sbras.nsc.ru/interval/Library/InteBooks/ SharyBook.pdf
[51 Шарый С.П. Оптимальное внешнее оценивание множеств решений интервальных систем уравнений. Часть 1 /7 Вычиел. технологии. 2002. Т. 7, № 6. С. 90 113.
Шарый С.П. Оптимальное внешнее оценивание множеств решений интервальных систем уравнений. Часть 2 /7 Там же. 2003. Т. 8, № 6. С. 84 109.
[6] SllARY S.P. A new class of algorithms for optimal solution of interval linear systems /7 Interval Comput. 1992. No. 2(4). P. 18 29.
[7] Ronx J. Systems of linear interval equations /7 Linear Algebra and its Appl. 1989. Vol. 126. P. 39 78.
[8] Beeck H. Uber die Struktur und Abschätzungen der Lösungsmenge von linearen Gleichungssystemen mit Intervallkoeffizienten /7 Computing. 1972. Vol. 10. P. 231 244.
[9] Nickel K. Die Überschätzung des Wertebereiches einer Funktion in der Intervallrechnung mit Anwendungen auf lineare Gleichungssysteme /7 Ibid. 1977. Vol. 18. P. 15 36.
[10] Oettli W. On the solution set of linear system with inaccurate coefficients /7 SIAM J. Numer. Anal. 1965. Vol. 2. P. 115 118.
[11] Задачи линейной оптимизации с неточными данными / М. Фидлер, И. Недома, Я. Рамик, И. Рои, К. Циммерманн. М. Ижевск: РХД, 2008.
[12] Krawczyk K.R., Neumaier A. An improved interval Newton operator /7 J. of Math. Analisis and Appl. 1986. Vol. 118. P. 194 201.
[13] Mayer. G. Epsilon inflation in verification algorithms /7 .J. of Comput. and Appl. Math. 1995. Vol. 60. P. 147 169.
[141 l^i Ml- S.M A note on epsilon-inflation /7 Reliable Comput. 1998. Vol. 4. P. 371 375.
[15] Neumaier A. New techniques for the analysis of linear interval equations /7 Linear Algebra and its Appl. 1984. Vol. 58. P. 273 325.
[16] Neumaier A. A simple derivation of Hansen Bliek Rolm Ning Kearfott enclosure for linear interval equations /7 Reliable Comput. 1999. Vol. 5. P. 131 136.
[17] Hargreaves G.I. Interval analysis in MATLAB /7 Manchester Center for Computational Mathematics, 2002. (http: //old. ict. nsc. ru/interval/Programming/INTLABtutor. pdf)
[18] Бауэр Ф.Л., Гооз Г. Информатика. М.: Мир, 1996.
[19] Панков П.С. Алгоритмы доказательства устойчивых утверждений и глобальной оптимизации в охраниченной области. Фрунзе, 1984 (Деп. в ВИНИТИ, № 5250-84Деп).
[20] Neumaier A. Interval methods for systems of equations. Cambridge Univ. Press, 1990.
[21] Toft O. Sequential and parallel solution of linear interval equations /7 Eksamensproject: NI-E-92-04. Numerisk Institute, Danmarks Tekniske Ho.jskole. Lyngbv, 1992.
[22] Gregory R.T. A collection of matrices for testing computational algorithms. N.-Y.: Wilcy-Intersci., 1978.
[23] Verification software in MATLAB / INTLAB. Электронный ресурс, доступный на http://www.cs.cas.cz/rohn/matlab
Поступила а редакцию 4 мая 2011 г., с доработки 16 сентября 2011 г.
