Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. Под ред. А.А.Абрамова (1986) (1095855), страница 49
Текст из файла (страница 49)
Если теперь положить Д=(1 — 21ч,и', )(1 — 2~'ггчг )... (1 — 2 1ч„,в'„г). (6,3ЛЗ) то (6.3.12) можно переписать как Д А =Я, Так как каждая матрица т 1 — 2и~;в~ ортогональна, то и их произведение Д также будет ортогональной матрицей (см. упражнение 6.34). Следовательно, Д ' =Д~, и (6.3.12) эквивалентно соотношению А=ДА, «6.3.14) которое и представляет собой ДЛ-разложение матрицы А. Подчеркнем, что это разложение может быть осуществлено. всегда, без каких-либо ограничений на матрицу А, Кроме того, указанный алгоритм разложения численно устойчив, и при его выполнении не может возникнуть никаких серьезных проблем, связанных с ошибками округления.
Мы снова без доказательства приводим теорему о ДЯ-алгоритме, аналогичную теореме 6.3.1, Т е, о р е м а 6.3.2, Пусть выполнены условия 2 и 3 теоремы 6.3.1. Тогда генерируемые ЯЯалгоритмом (6.3.9) матрицы А» сходятся к верхней треугольной матрице вида (6.3.8). При этом лежащие ниже главной диагонали элементы а,(~~1 матриц А ь сходятся к нулю с -линейной скоростью, определяемой отйошснием собственных значений матрицы А, а именно а~ ~ =0~ 1 ~Л,!" (6.3.15) б ~ ~1 ~lс Эта теорема является теоретическим обоснованием ( Ж-алгоритма. Однако как в формулировке теоремы, так и в схеме алгоритма имеется ряд ограничений, которые снижают практическую ценность метода. Мы сейчас рассмотрим эти ограничения и обсудим способы, которые позволяют их обойти, по крайней мере частично.
я -~, г)1 Обратимся сначала к условию 2 теоремы 6.3.1, которое требует, чтобы все собственные значения матрицы А бьши различны по абсолютной величине. В случае вещественной матрицы это исключает не только кратные собственные значения, но и пары комплексно-сопряженных собственных значений. Однако на самом деле ситуация не настолько плоха. Если матрица А вещественна, то все множители Д~ и Я~ также будут вещественны, и матрицы Аь, разумеется, не смогут сходиться к треугольной матрице с невещественными собственными значениями. Тем не менее оказывается (и это самое большее, на что мы вообще могли надеяться), что Аь "сходятся"* к почти треугольной форме, структура которой иллюстрируется мат- рицей Ха е ~9 Мы здесь предположили, что матрица имеет три вещественных собственных значения Хз, Ха, Х9 с различными абсолютными величинами и три пары комплексно-сопряженных собственных значений также с различными абсолютными величинами.
Эти пары собственных значений определяются тремя 2 Х 2-матрицами, которые показаны в виде блоков на главной диагонали. На самом деле элементы соответствующих 2 Х 2-матриц последовательности Аь могут не сходиться, однако собственные значения этих матриц будут сходиться к собственным значениям А. Таким образом, в рассматриваемом случае вычисление комплексных собственных значений вещественных матриц не представляет проблемы.
ДЛ-алгоритм в том виде, как он описан выше, является очень неэффективным. Мы сейчас рассмотрим две важных модификации, которые необходимо внести в алгоритм, чтобы его можно было использовать на практике. Первая проблема связана с тем, что разложение (6.3.14) требует 0(п') операций, так что каждый шаг процесса выполняется слишком медленно. Мы можем обойти эту трудность, если предварительно преобразуем матрицу А к форме, которая позволяет проводить это разложение за меньшее время.
На самом деле нам бы хотелось сначала привести матрицу к трехдиагональной форме, как это делалось в предыдущем разделе, а затем к полученной трехдиагональной матрице применить описаннный ДЯ-алгоритм. К сожалению, хотя такое приведение к трехдиагональной форме и возможно, в случае несимметричной матрицы оно может оказаться численно неустойчивым. Самое большее, чего мы можем достичь с гарантией численной устойчивости, это привести матрицу к виду л11 я!3 (6.3.16) А= 0 1 яп д соз Π— созд япд (6.3.17) где первый синус стоит в позиции с индексами ( 1. 1) . Умножая теперь матрицу Хессенберга (6.3.16) слева на матрицу вращения Р,, у которой синусы и косинусы расположены в левом верхнем, углу, получим ялд совд — сов д яп О 1 а1„ а,, авн ни,п а1н яп О +аз„сов д -а1н совд+из„зшд а,1яп О+а,,совд — а,, совд+ аз, яп О аз авв Такие матрицы, у которых ниже главной диагонали имеется только одна ненулевая диагональ, непосредственно примыкающая к главной, называются матрицами Хессенберга.
Исходная матрица А может быть эффективно приведена к форме Хессенберга с помощью описанных в предыдущем разделе преобразований Хаусхолдера. Действительно, если определить векторы и, (1 = 1,..., п — 2) точно так же, как в разделе 6.2, и выполнить преобразования подобия (6.2.1) с ортогональными матрицами Р~ =1 — 2иг;и т то в результате мы получим матрицу вида (6.3.16).
В упражнении 6.3.6 предлагается доказать это утверждение. Теперь будем считать, что матрица А приведена к форме Хессенберга, и 'рассмотрим вопрос о применении к ней ДА-алгоритма. Прежде всего отметим, что в этом случае разложение (6.3.14) осуществляется особенно просто. Это можно было бы сделать, как н раньше, с помощью преобразований Хаусхолдера, но несколько лучше воспользоваться' операциями, известными как преобразования Гивенса, или плоские вращения.
Эти операции задаются ортогональными матрицами вида Если выбрать угол 0 из условия — а,, созО+аз, яп8=0, или О=агстй(а,,/а~,), то элемент - (2,1) произведения обратится в нуль. Таким обра- зом, в результате умножения на последовательность матриц плоских вращений Р; с соответствуюшими значениями углов все поддиагональные элементы матрицы Хессенберга последовательно обратятся в нуль. Итак, мы приходим к разложению А =1„1Л, где Д =Р„, ...Р,, и, следова- тельно, 0=Р, ...Р~ ~. Это разложение может быть получено за 0(л') операций, что следует считать достаточно эффективным. Приведение матрицы А к форме Хессенберга не имело бы смысла, если бы эту процедуру пришлось выполнять после каждого шага алгоритма.
К счатью, это не так. Действительно, отличные от нуля внедиагональные элементы матриц Р, находятся только в позициях (/ + 1, /) и (/,/ + 1). Отсюда легко следует (упражнение 6.3.7), что матрица Д сама является матрицей Хессенберга. Далее, так как Я вЂ” верхняя треугольная матрица, произведение ЯД также является матрицей Хессенберга. Таким образом, если исходная матрица приведена к форме Хессенберга, то и все матрицы А~„генерируемые Дй-алгоритмом, будут автоматически сохранять эту форму. Даже при первоначальном приведении к форме Хессенберга и вытекаю- щей отсюда экономии при построении матриц А ~ алгоритм может оказать- ся неэффективным из-за медленной сходимости к нулю элементов, находя- щихся под диагональю. Скорость сходимости определяется выражением (б.3.15), нз которого видно, что если два собственных значения, скажем, Л~ и Л~+~, близки, то внедиагональный элемент с индексами (/ + 1,/) бу- дет стремиться к нулю очень медленно.
Попытаемся ослабить эту проблему сходимости следуюшлм образам. л Предположим, что Л„является хорошим приближением к Л„, и рассмотл л рим матрицу А = А — Л„/, собственные значения которой равны Л, — Л„,... л л ..., ˄— Л„. Если теперь применить ДА-алгоритм к А, то внедиагональные элементы последней строки матриц А ~ будут стремиться к нулю как стел л пени отношений (˄— Л„)/(Л; — Л„), а не как степени Л„/Л; (/ = 1,...
..., и — 1) . Наиболее медленно из ннх сходится элемент в позиции(п, п-1), л причем для матрицы А скорость сходимости определяется отношением л л (˄— Л„) /(Л„, — Л„), а для исходной матрицы А — отношением Л„/Л„ л Если, например, Л„=0,99, Л„~ =1,1 и Л„= 1,0, то Л„/Л„~ =0,9, а ! ˄— — Л„! /! Л„, — Л„! = 0,1, так что элемент (и,п — 1) матрицы со сдвину- тыми собственными значениями сходится к нулю примерно в 20 раз быст- рее. Хорошее приближение Л„, которое следовало бы использовать в качестве параметра сдвига, нам обычно неизвестно. Однако в ходе ДЯ-алгоритма элементы а„„матриц Аь будут сходиться к Л„, так что мы можем их И1 использовать как параметры сдвигов, т.е., сделав Й шагов, мы можем на шаге /с+ 1 выполнить ДЯ-разложение матрицы А» =Ах — а„„Е В то же л 1У,.1 времяможно выполнять преобразование сдвига на каждом шаге Дй-алго- ритма, используяв качестве параметра элемент (л,л) текущей матрицы. ?06 Так как при каждом сдвиге собственные значения исходной матрицы изменяются на величину параметра сдвига, то необходимо следить занакоплением этих величин.
Фактически, именно сумма сдвигов сходится к собственному значению Х„. Критерием сходимости служит достаточная малость элементов последней строки. Когда эта малость будет достигнута, можно отбросить последнюю строку и последний столбец матрицы н перейти к определению собственного значения Х„,, исходя из полученной подматрицы размера (н — 1) Х (и — 1) . Отметим, что собственные значения этой подматрицы, а следовательно, и исходной матрицы, были изменены на суммарную величину всех сдвигов (которая служила приближением к Х„), так что после вычисления собственных значений подматрицы к ним следует добавить эту величину.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
