Дж. Деммель - Вычислительная линейная алгебра, страница 71

312 Глава б. Итерационные методы для линейных систем Алгоритм 6.8. Симметпричный метод БОВ (ББОВ/: 1. Выполнить шаг метода БОК(ы), вычисляя компоненты вектора х в обычном порядке возрастания х; а, х;д, ..., х; „. 2. Выполнить шаг метода БОВ(ы), вычисляя компоненты вектора х в обратном порядке х, „, хса г,, х;,г. Мы представим этот алгоритм в виде хгьг = Е х; + с и покажем, что матрица Е имеет вещественные собственные значения, а потому применение чебышевского ускорения становится возможным. Предположим, что А — симметричная матрица (как это имеет место в модельной задаче), и представим А, как в уравнении (6.19): А = 0 — Х вЂ” Р = Ю(1 — 1 — 1/).

Поскольку А = Ат, имеем 1/ = Ьт. С помощью равенства (6.21) перепишем два полушага метода ББОВ. следующим образом: 1. Хг ьа/г — (1 — юЬ) ((1 — ш)1 + О/Г/)хв + са/2 = Ьмх1 + са/г~ 2. хгьг = (1 — юс/) ((1 — ы)1 + ы1)хы.~/г + сг = (/ух~+~/г + сг. Исключая отсюда х;+,/г, получим хеьг = Е х, + с, где Е = 1/„Ь 1 + (ы 2)г(1 ,ц[/)-~(1 ~,Ь)-~ + (ы 2)(1 ы(/)-~ + (ы — 2)(1 — ыУ) ~(1 — ыЬ) г(1 — ы(/). Мы утверждаем, что собственные значения матрицы Е„вещественны. Действительно, она имеет те же собственные значения, что и подобная ей матрица (1 — ш(/)Е (1 — ы(/) = 1+ (2 — ы)г(1 — ы1) '(1 — ы(/) г + (ы — 2)(1 — ы(/) + (ы — 2)(1 — ыЬ) — 1 + (2 ,)г(1 ,1,)-~(1 ,„1,т)-~ + (, 2)(1 1,т) -а + (ы — 2)(1 — огЬ) Очевидно, что эта матрица симметрична и потому должна иметь вещественные собственные значения.

Пример 6.12. Применим метод ББОВ(ы) с чебышевским ускорением к модельной задаче. Нужно выбрать значение ы и оценить спектральный радиус р = р(Е ). Оптимальное значение ы, минимизирующее р, не известно, однако Янг показал (см. [267, 137]), что вполне удовлетворителен выбор —,-.Хел ~- ..и,.---- бышевского ускорения погрешность на шаге т умножится на число р га — (-,' — -) < 2/(1+т~Я). Поэтому потребуются т = 0(а/'/г) = 0(п'/~) итераций, чтобы уменьшить погрешность в фиксированное число раз, большее единицы. Поскольку каждая итерация имеет ту же стоимость 0(п), что и итерация метода БОВ(ы), общая стоимость процесса составляет 0(па/4). Это объясняет элемент табл. 6.1, относящийся к методу ББОВ с чебышевским ускорением. В отличие от сказанного выше, после т шагов метода БОВ(ы,рг) погрешность умножилась бы лишь на (1 — —;,) .

К примеру, положим А/ = 1000. Тогда, 6.6. Методы крыловского подпрострвпствв чтобы уменьшить погрешность вдвое, методу ВОВ.(ы, ~) требуются 220 итераций, тогда как для %0В(ы,р~) с чебышевским ускорением нужны лишь гп = 17 итераций. О 6.6. Методы крыловского подпрострапства Эти методы используются и для решения системы Ах = Ь, и для вычисления собственных значений матрицы А. В них предполагается, что А доступна лишь посредством ечерного ящика» в виде подпрограммы, вычисляющей по заданному вектору в произведение у = Ав (а также, возможно, произведение у = А~в, если матрица А несимметрична). Другими словами, прямой доступ к элементам матрицы и их изменение не используются.

Такое предположение разумно по нескольким причинам. Во-первых, умножение на вектор— это наиболее дешевая нетривиальная операция, которую можно проделать с (разреженной) матрицей. Если в А имеется т ненулевых элементов, то для матрично-векторного умножения нужны тп скалярных умножений и (самое большее) т сложений. Во-вторых, А может не быть представлена явно в вцде матрицы, а доступна лишь через подпрограмму для вычисления произведений Ах. Пример 6.13. Предположим, что имеется физическое устройство, поведение которого моделируется программой, вычисляющей по вектору х входных параметров вектор у выходных параметров, описывающих поведение устройства. Выходной вектор у может быть как угодно сложной функцией у = 7(х), возможно, требующей даже решения нелинейных дифференциальных уравнений. Например, параметры, составляющие вектор х, могут описывать форму крыла, а 7'(х) может быть лобовым сопротивлением на крыле, вычисляемым путем решения уравнений Навье — Стокса для потока через крыло.

Типовой конструкторской задачей является выбор входа х, оптимизирующего поведение устройства 7(х); для определенности, предположим, что оптимизация означает наибольшее возможное уменьшение 7(х). Наша задача в этом случае сведется к тому, чтобы насколько возможно удовлетворить уравнение 7(х) = О. Для целей иллюстрации, предположим, что х и у — векторы одинаковой размерности.

Тогда очевидным выбором для приближенного решения уравнения является метод Ньютона, что приводит к итерациям х(~"') = х<~> — (~7(х<~))) '7(х~ ~), где '7,7(х~~>) — матрица Якоби функции 7" в точке х( 1. Это соотношение можно переписать как задачу решения линейной системы (т7(х~™~))б< ~ = 7(х(™>) относительно б< ~ с последующим вычислением вектора х< +О = х1 > — б< ).

Как же, однако, следует решать эту линейную систему с матрицей коэффициентов '7,7(х< ~), если даже вычислением(хпв>) представляет трудность? Оказывается, что для произвольного вектора г можно вычислить матрично-векторное произведение (~77'(х)) в, поэтому к решению линейной системы можно применить методы крыловского подпространства. Используя для вычисления (177'(х)) в разностные оглношенил или формулу Тейлора, можно видеть, что (7(х+йв) — 7(х))/й (~77(х))% Таким образом, чтобы вычислить (177"(х)) в, нужно дважды обратиться к подпрограмме, вычисляющей 7( ), один раз с аргументом х, а другой — с аргументом х + Ьв.

Иногда, однако, бывает трудно выбрать значение и, дающее 314 Глава б. Итерационные методы для линейных систем хорошую аппроксимацию производной (выбор слишком малого Ь приводит к потере точности из-за округлений). Другой способ вычисления (~77" (х)) х состоит в реальном дифференцировании функции 7". Если 7" — достаточно простая функция, то это можно проделать на бумаге. Для сложных функций 7" существуют инструментальные средства, позволяющие по (почти) любой подпрограмме, вычисляющей 7(х), автоматически породить другую подпрограмму, вычисляющую ('77(х)) е (см.

[29]). Можно также, хоть это и менее эффективно, использовать специальные возможности языков С++ и фортран 90. О Существует множество различных методов крыловского подпространства. Некоторые из них пригодны для несимметричных матриц, другие требуют от матрицы симметрии или положительной определенности.

Некоторые методы для несимметричных матриц предполагают, что могут вычисляться не только произведения типа Ах, но и произведения типа Атх; в зависимости от способа представления матрицы А, такие произведения могут быть или не быть доступны (см. пример 6.13). Наиболее эффективным и лучше всего понимаемым методом этого класса является метод сопряженных градиентов (метод СО). Он применим только к симметричным положительно определенным матрицам, в частности он приложим к модельной задаче. В этом разделе мы сосредоточимся на методе СС. Если матрица не является симметричной положительно определенной, то выбор наилучшего метода из многих существующих может представлять определенные трудности.

В разд. 6.6.6 дан краткий обзор методов, имеющихся помимо метода СС, а также рекомендации по поводу того, какой метод в какой ситуации следует использовать. Отсылаем читателя также к более полному сетевому справочнику в НЕТЫВ/1ешр1аФеэ, включающему в себя книгу [24] н реализации методов на языках МаНаЬ, фортран и С++. Обзоры современных исследований по методам крыловского подпространства можно найти в [15, 107, 136, 214!.

В гл. 7 мы обсудим методы крыловского подпространства в применении к задаче вычисления собственных значений. 6.6.1. Извлечение информации о матрице А из матрично-векторных произведений Пусть даны вектор Ь и подпрограмма, вычисляющая произведения вида А аз что можно узнать с их помощью о матрице А7 Наиболее очевидным решением является вычисление последовательности матрично-векторных произведений дг = Ь, уз = Ауы уз = Ауг = А уы ..., у„= Ау„~ = А" 'уы где п — порядок матрицы А. Положим К = [уы ую..., д„]. Тогда можно написать А К = [Ауы,Ау„мАуо] = [уэ у„,А"уг] (629) Заметим, что первые и — 1 столбцов матрицы А К совпадают с последними п — 1 столбцами в К (при нумерации, сдвинутой на единицу). Предположим пока, что матрица К невырожденна; тогда можно вычислить вектор с = — К зА"уь Теперь А.

К = К. [ею ее,, е„, — с] н К С, б.б. Методы крыловского подпространства где е; — столбец 1 в единичной матрице. Отсюда О О .. О 1 О О О 1 К 'АК=С= 1 — с„ Обратим внимание, что С вЂ” верхняя хессенбергова матрица. В действительности, это сопровоокдающая матрица (см. равд.

4.5.3) многочлена р(х) х" + 2 ',ы, с,х' ~, ЯвлЯющегосЯ ее хаРактеРистическим многочленом. Таким образом, используя лишь матрично-векторные произведения, мы привели А к очень простой форме, позволяющей, в принципе, найти собственные значения как нули многочлена р(х). Однако эта простая форма бесполезна для практики в силу следующих причин: 1. Определение вектора с требует п — 1 матрично-векторных умножений с участием матрицы А, а затем решения линейной системы с матрицей К. Даже если А в разреженная матрица, К, скорей всего, будет плотной, поэтому нет никаких оснований ожидать, что линейная система с матрицей К решается сколько-нибудь проще, чем исходная система Ах = 5.

2. С большой вероятностью, матрица К будет очень плохо обусловленной, а потому с будет вычислен очень неточно. Причина в том, что для вычисления столбцов рл матрицы К наш алгоритм выполняет степенной метод (алгоритм 4.1), поэтому векторы р, сходятся к собственному вектору, отвечающему наибольшему собственному значению матрицы А. Таким образом, столбцы в К все более и более приближаются к параллельным. Мы преодолеем эти проблемы следующим образом: заменим К ортогональной матрицей Я, такой, что при любом к линейные оболочки первых к столбцов в Х и Я являются одним н тем же подпространством. Оно называется крыловским подпространством. В отличие от К, матрица Ч хорошо обусловлена и легко обратима.