ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 11, с. 1811-1818
УДК 519.614.2
МОДИФИКАЦИЯ НЕКОТОРЫХ МЕТОДОВ ТИПА СОПРЯЖЕННЫХ НАПРАВЛЕНИЙ ДЛЯ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ1)
© 2007 г. Л. Ф. Юхно
(125047 Москва, Миусская пл., 4а, ИММ РАН) e-mail: yukhno@imamod.ru Поступила в редакцию 26.03.2007 г.
Предлагается и исследуется модификация некоторых известных методов типа сопряженных направлений. Полученные модифицированные методы более устойчивы относительно накопления погрешностей арифметических операций. Эти методы применимы для решения плохо обусловленных систем линейных алгебраических уравнений, возникающих, в частности, при аппроксимации некорректных задач. Приводятся результаты численных экспериментов, иллюстрирующие преимущество предложенных модификаций. Библ. 4. Табл. 2.
Ключевые слова: система линейных алгебраических уравнений, метод сопряженных направлений, численная устойчивость.
1. ВВЕДЕНИЕ
Общим свойством методов типа сопряженных направлений для решения систем линейных алгебраических уравнений (СЛАУ) является разложение искомого решения по совокупности векторов, ортогональных по некоторому скалярному произведению. Ненулевых таких векторов может быть не более размерности пространства, которому принадлежат все эти векторы. Поэтому каждый такой метод при точной реализации завершается за конечное число шагов, т.е. является конечным методом.
Важным этапом реализации этих методов является ортогонализация выбираемой совокупности векторов по соответствующему скалярному произведению. Как известно (см. [1]), при проведении процесса ортогонализации, особенно в задачах с плохо обусловленной матрицей, может произойти значительная потеря точности. В результате процесс построения "ортогональной" системы векторов для некоторых таких методов становится бесконечным, т.е. эти методы можно рассматривать как чисто итерационные (см. [2]). Однако поскольку такое итерационное продолжение полностью определяется погрешностями вычислений, то иногда в процессе итераций может начаться "разбалтывание", т.е. увеличение значения функционала ошибки, соответствующего данному методу, при переходе к последующим приближениям, в то время как это значение теоретически не должно возрастать. Это подтверждается численными экспериментами для задач с плохо обусловленными матрицами.
В работе предлагается модификация некоторых методов такого типа, а именно рассматриваются методы, приведенные в [1] в главе "Метод минимальных итераций и другие методы, основанные на идее ортогонализации". Модифицированные методы при точных вычислениях полностью совпадают с исходными, а при неточных - по той же причине становятся итерационными. Однако итерации строятся таким образом, что функционал ошибки на каждом шаге не возрастает. Поэтому хотя предлагаемые модификации также подвержены вредному влиянию погрешностей округления при реализации процесса ортогонализации, указанного разбалтывания в этих методах быть не должно, т.е. их устойчивость гарантирована.
В работе приводятся результаты численных экспериментов, иллюстрирующие сравнительные характеристики исходного и модифицированного вариантов некоторых из этих методов на модельных примерах.
Настоящая работа является продолжением работы [3], где предложена модификация метода сопряженных градиентов после первой и второй трансформации Гаусса для решения СЛАУ, а также метода решения задачи исключения, предложенного в [4].
^Работа выполнена при финансовой поддержке РФФИ (код проекта 06-01-00749).
1811
1812
ЮХНО
2. ПОСТАНОВКА ЗАДАЧИ И СУТЬ ПРЕДЛАГАЕМОЙ МОДИФИКАЦИИ Будем рассматривать СЛАУ
Ах = Ь, (2.1)
где А - в общем случае матрица порядка т х п, вектор Ь е С™ задан, х е С" искомый. Для пары векторов как из С", так и из С™ будем использовать стандартное скалярное произведение
( и, V) = V* и.
Предполагается, что уравнение (2.1) имеет решение, может быть неединственное. Обозначим через х* то решение, которое получается рассматриваемым методом. Для этого метода фиксируется некоторая эрмитова положительно-определенная п х п-матрица Я. При х е С" функция ^(х) = (Я(х - х*), х - х*) есть функционал ошибки. Как указано в [1], при точной реализации таких
методов значения соответствующих им функционалов ошибки не возрастают на каждом шаге процесса.
Общим для всех рассматриваемых методов является то, что последующее приближение х может быть получено из предыдущего х по формуле
х = х + ag.
Здесь а - число, а g - некоторый вектор, соответствующий данному методу. При этом соответствующее значение а выбирается из условия, что х) = Е(х + ag) как функция а достигает при этом значении минимума.
Идея модификации состоит в том, что во всех рассматриваемых далее методах вектор Rg представляется в виде
Rg = А* д, (2.2)
где д - вектор, получаемый непосредственно при реализации метода или несложными дополнительными вычислениями. В этом случае (см. вывод в [3]) значение а, дающее минимум функционала ошибки Е(х + ag), определяется формулой
а = (Ь - Ах, д) = (Ь, д) - (х, Rg)
а ( яё, g) ( g) ' (2.3)
которая представляет суть модификации. В этой формуле в отличие от исходных вариантов методов учитывается отклонение численно получаемых векторов от ортогональности.
При точной реализации вычислений модифицированная формула (2.3) эквивалентна исходным вариантам этой формулы. Это легко может быть проверено непосредственно для каждого метода.
При неточной реализации вычисление значения а по формуле (2.3) не влечет существенных погрешностей. Это гарантирует, что х) ^(х) независимо от того, как было получено g, т.е. функционал ошибки не возрастает на каждом шаге процесса. В силу этого модифицированный метод позволяет продолжать итерации для к > п с целью получения большей точности. Конечно, накопление погрешностей арифметических операций имеет место и в этом случае, однако оно не является "катастрофическим" для метода, так что разбалтывание можно считать исключенным.
Таким образом, использование для вычисления а формулы (2.3) вместо исходной при условии выполнения (2.2) и составляет предлагаемую модификацию рассматриваемых методов. Далее рассматриваются конкретные методы и соответствующие им модификации. В настоящей работе мы ограничиваемся случаем т = п, detA Ф 0, что гарантирует существование и единственность решения х* системы (2.1).
3. МЕТОД ОРТОГОНАЛЬНЫХ ВЕКТОРОВ
Метод применяется для случая, когда А - эрмитова положительно-определенная матрица. В качестве Я берется матрица А. В этом методе система ортогональных и нормированных векторов исходного базиса въ ..., еп подвергается А-ортогонализации. В результате получаются век-
торы ...,gn, где
= е1, gi = ег - У г1 gl- ••• - У г, г-^г -1, (Аег, gj) ~ - ■ 1 т -1
У г] = ГИТ^Тл, г = 2 П' Л = 1 г -1-Параллельно вычисляется последовательность векторов хг, г = 0, 1, ..., п, по формулам
х0 = 0, хг = хг-1 + а^ при г = 1, 2, ..., п,
где
а = (б, gi) / (Agi, gi). (3.1)
В результате получается х* = хп.
Как отмечено в [1], А г = Ае г _ это 1-й столбец матрицы А, а также
(Agi, gi) = (Аег, gi) = (Аг, gi). (3.2)
Так как в рассматриваемом методе Л = А = А*, то для выполнения равенства (2.2) возьмем q = g, т.е. в соответствии с рекомендациями разд. 2 используем вместо формулы (3.1) формулу
а = (Ь - Ахг - 1, g1 ) = (Ь, gг ) - ( хг - 1, Аg1 ) (3 3)
г (Agг, gi ) (Agг, gi ) . (3.3)
Формулы (3.1) и (3.3) эквивалентны при точной реализации вычислений. Действительно, так как хг _ 1 - линейная комбинация векторов g1, ..., g1 _ 1, а каждый из них А-ортогонален вектору g1, то дополнительное слагаемое (х г_ 1, Ag1) в формуле (3.3) равно нулю. При наличии погрешностей это слагаемое ввиду нарушения ортогональности векторов нулю не равно и указанные формулы могут дать различающиеся результаты. При этом использование формулы (3.3) гарантирует, что
(А ( х * - хг), х * - хг ) ^ (А ( х * - хг -1 ), х * - хг -1 ), (3.4)
а использование формулы (3.1) такой гарантии не дает.
Отметим также, что равенство (3.2) тоже предполагает точную реализацию процесса А-орто-гонализации. Учитывая упомянутую в разд. 2 неприятность, при решении плохо обусловленных задач для вычисления (Agi■, gi) лучше этой формулой не пользоваться.
Таким образом, предлагаемая модификация рассматриваемого метода включает использование формулы (3.3) вместо (3.1), а также неиспользование упрощающей формулы (3.2). В результате гарантируется выполнение соотношения (3.4) при вычислениях с погрешностями.
4. МЕТОД А-МИНИМАЛЬНЫХ ИТЕРАЦИИ
В этом методе А _ эрмитова положительно-определенная матрица, Я = А. Для произвольно взятого вектора 5 рассматривается последовательность векторов 5, А.^, А25, ..., Ап_ ^ которая предполагается линейно независимой. Эта последовательность векторов подвергается процессу А-ортогонализации. Здесь последовательность векторов g1, ..., gn строится по следующим трехчленным рекуррентным формулам:
gl = 5, g2 = Agl _ ч&ъ gi + 1 = Agi _ у& _ _ 1 при г = 1, 2, ..., п _ 1,
где
(Agг, Agi) 8 (Agг, gi)
Уг = / л „ „ - , 8г =
( Agг, &)' г ( Agг-1, gг _1 )'
Последовательность х0, ..., хп строится по формулам
х0 0, х1 +1 хг
1814
ЮХНО
где
a- = (М) • ' -1-2- "<4Л)
Как и в разд. 3, возьмем q = g, т.е. заменим формулу <4.1) формулой
a = (b - Axr, gr) = (b, g,) - ( x - Agi) <42)
a ( Agi, g- ) ( Agi, g-) • <4.2)
Аналогично предыдущему можно показать эквивалентность формул <4.1) и <4.2) при их точной реализации, а также преимущество использования модифицированной формулы <4.2) перед <4.1) для решения плохо обусловленных задач, когда влияние погрешностей округления значительно.
5. МЕТОД СОПРЯЖЕННЫХ ГРАДИЕНТОВ
В этом методе A - эрмитова положительно-определенная матрица, R = A. Расчетные формулы метода таковы:
Хо = 0, Го = b • g| = r0; далее при i = 1, 2, ..., n вычисляются
(r--|, ri_|) (ri_|, gi) (Го, gi)
ai =
( Agi, gi ) ( Agi, gi ) ( Agu gi У х- = х- _| + aigi ;
при i = 1, 2, ..., n - 1 вычисляются также
Гi = Гi _|- ai Agi,
e = ( r i, ri ) = _( г-, Agi)
Pi ( г- _ |, г- _ | ) ( Agi, gi ) •
gi + i = ri + P-gi •
В [1] для непосредственного вычисления a,
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.