1626435587-55f52a4de97976f3c6215fa7c103f544 (Смирнов 2017 - Основы вычислительной физики ч2), страница 5

Пренебрегая квадратичными членами(+1) u(+1) , получаем линейную систему на ( + 1)-й итерации:(︀)︀(︀)︀ − () u(+1) − (+1) u() = − − () u() .Поскольку собственный вектор определён с точностью до множителя,(+1)можно положить = 0, в результате чего получим систему из (+1)(+1)уравнений на неизвестных величин 1, . . . , −1 , (+1) :22⎛(+1)1⎜..⎜.

() · ⎜⎜⎝ (+1)−1(+1)⎛11 − ()⎜..⎜.≡⎜⎜⎝ −1,1,1⎞⎟⎟(︀)︀⎟ = − − () u() ,⎟⎠12...−1,2,2...где () ≡1,−1.... . . −1,−1 − ()...,−1()−1...()−−1()−⎞(40)⎟⎟⎟⎟⎠.Решение исходной спектральной задачи получается в результатеитерационного процесса u(+1) = u() + u(+1) , (+1) = () + (+1) .1.7.5. Метод вращений ЯкобиОднократное применение большинства рассмотренных выше методов позволяет найти лишь одно собственное значение и соответствующий ему собственный вектор матрицы. Для полноты картины рассмотримдля одновременного нахождениясобственных значений и векторов симметричных вещественных матриц6 .

Идея метода Якоби лежит в основе большинства используемых внастоящее время наиболее эффективных алгоритмов: матрица приводится к диагональному виду Λ путём последовательности преобразований подобия.Действительно, из курса алгебры известно, что если — симметричная матрица, то существует ортогональная матрица ( = −1 ),такая что −1 = Λ, где Λ — диагольная матрица, состоящая изсобственных значений матрицы .

В методе Якоби задача нахождения преобразования подобия решается методом последовательныхприближений: строится последовательность матриц (0) = , (1) ,(2) , . . ., так что () → Λ. Для построения следующего приближения (+1) по заданной матрице () необходимо вначале найти максимальный по модулю элемент, лежащий выше главной диагонали, иметод вращений Якобивсех6 Помимо того, что симметричные (эрмитовы) матрицы играют важную роль вомногих областях физики, эрмитовость матрицы является достаточным условиемустойчивости собственных векторов по матричным элементам, в то время как спектральная задача для несимметричных матриц может быть плохо обусловленной.Красивый пример высокой чувствительности несимметричной матрицы Уилкинсона к малому возмущению матричных элементов приведён в монографии [2, с. 161].Более подробный анализ вопроса можно найти в книге [A3, с. 93].23запомнить его индексы , : , :⃒ () ⃒⃒ () ⃒⃒⃒⃒⃒ = max .(41)<Основная идея метода заключается в том, чтобы с помощью поворотав плоскости ( , )(42)(+1) = ( ) () ( )занулить максимальный по модулю недиагональный элемент ( , ) наследующем, ( + 1)-м, шаге итерационного процесса, приблизив темсамым матрицу (+1) к искомой диагональной матрице Λ.

На следующей итерации данный элемент вновь будет отличен от нуля, однаконесложно увидеть, что после каждого поворота (42) сумма квадратоввнедиагональных элементов матрицы будет уменьшаться на величину 22 , а сумма квадратов диагональных элементов, напротив, возрастёт на ту же величину: сумма квадратов всех матричных элементов есть скаляр и не изменяется при поворотах системы координат.Отсюда очевидно, что итерационный процесс, построенный на занулении максимальных по модулю недиагональных элементов посредствомплоских вращений, будет сходиться к искомой диагональной матрицеdiag{1 , . . . , }.

Результирующая матрица поворота, равная произведению матриц плоских вращений на каждом шаге, будет содержатькоординаты собственных векторов. Рассмотрим детали алгоритма иформулы для эффективного преобразования матрицы .Поворот в плоскости ( , ) описывается ортогональной матрицей( ) , для элементов которой можем записать: = 1 при ̸= , ̸= ; = 0 при ̸= , ̸= , ̸= , ̸= ; = = ;− = = ,(43)где для краткости использованы обозначения ≡ cos ,(44) ≡ sin .Подставляя (43) в (42), получим для элемента ( , ) матрицы (+1)()()(+1):() = (2 − 2 ) + · · ( − ).Выбрав так, чтобы()tg(2 ) =2 () () −24,| | ≤,4(45)(+1)получим = 0, что приблизит матрицу (+1) к диагональной.Для выполнения (45) можно вычислить значения ≡ cos и ≡ sin по следующим формулам:√︃ (︂√︃ (︂)︂)︂11||||1+, = sign( )1−,(46)=22где = − ,=√︁ 2 + 42 .(47)Вычислив значения и по формулам (46) и (47) и подставивих в (43) и далее в (42), мы можем найти следующее приближение(+1) , сделав шаг итерационного процесса.

Заметим, однако, что непосредственное вычисление произведения матриц в формуле (42) крайненеэффективно ввиду большого количества нулей и единиц в ( ) (43).В самом деле, при повороте в плоскости ( , ) из 2 элементов матрицы () изменятся лишь 4 − 4 чисел в двух строках и двух столбцахс номерами и . В этой связи для выполнения плоского поворота(42) матрицы () следует отказаться от использования общего видаматричного умножения в пользу специальных формул (48–50), которые могут быть легко получены из (42) исключением суммированиянулевых элементов и умножений на единицу в многомерных матрицахплоского поворота. В соответствии с (48), бо́льшая часть7 коэффициентов матрицы () остаётся без изменений:(+1)()= при ̸= , ̸= , ̸= , ̸= .(48)(Греческими буквами , обозначены индексы, пробегающие в циклезначения от 1 до за исключением и ; в противоположность этому,индексы и определяются выражением (41) и являются фиксированными в пределах одной итерации.) Следующие две формулы определяют преобразование большинства коэффициентов в двух строках синдексами и :(+1)(+1)(+1)= =(+1) ()()при ̸= , ̸= ;() ,при ̸= , ̸= .= + ,=()−(49)Наконец, последние три формулы определяют преобразование коэф7В предположении большого размера матриц ≫ 1.25фициентов на пересечении строк и столбцов с номерами , :(+1) (+1)()()()()()()= 2 + 2 + 2 ; = 2 − 2 + 2 ;(︀ ()(+1)(+1)()() )︀ = = (2 − 2 ) + − .(50)Для удобства использования перечислим основные шаги методавращений Якоби:1) поиск максимального по модулю наддиагонального элемента вматрице () — определение индексов , (41);2) проверка условия малости | , |: если наибольший по модулюнедиагональный элемент матрицы () мал, завершаем итерационный процесс;3) вычисление = cos , = sin по формулам (46, 47);4) переход от матрицы () к (+1) по формулам (48–50);5) вычисление матрицы перехода = (0 0 ) (1 1 ) · .

. . · ( ) , столбцы которой при → ∞ дают собственные векторы матрицы .1.7.6. Сравнение методовСтепенной метод исключительно прост в реализации, однако достаточно эффективен лишь для ленточных матриц при |1 /2 | ≫ 1. В случае близких по модулю собственных значений скорость сходимости степенного метода падает обратно пропорционально разности 1 − |2 /1 |.Однако даже если собственные значения заметно различаются, другиерассмотренные выше методы, как правило, демонстрируют более высокую эффективность, что ограничивает практическую применимостьстепенного метода.Методы линеаризации и обратных итераций со сдвигом обеспечивают достаточно быструю сходимость (квадратичную вблизи корня)и особенно эффективны в случае ленточных матриц.

Дополнительнымпреимуществом по сравнению со степенным методом является возможность нахождениясобственного значения, для которого известно начальное приближение. Как и в случае трансцендентных уравнений, сходимость метода Ньютона не гарантируется: при неудачно выбранном начальном приближении итерационный процесс может бытьрасходящимся или сходиться не к искомому собственному значению, ак другому.Метод вращений Якоби годится лишь для симметричных вещественных матриц. Хотя он уступает в эффективности QR-алгоритму,любого26авторы монографии [4, с.

571] рекомендуют его к использованию приработе с не слишком большими матрицами ввиду высокой точности ипростоты программирования.Как и при рассмотрении других тем данного курса, в этой главе мыпознакомились лишь с узким кругом базовых методов, наиболее простых в понимании и реализации в программном коде. Вместе с тем данный раздел численных методов имеет и свою специфику. Спектральныезадачи, возникающие на практике, зачастую сопряжены с обработкойогромных массивов информации — матриц большого размера, возникающих при дискретизации непрерывных функций и операторов.

Вследствие этого вопросы эффективности используемых численных методов и их программных реализаций становятся критически важными,что заставляет задуматься об использовании программных решенийиз готовых библиотек. Авторы монографии [4, с. 567] пишут по этомуповоду: «Возможно, вы уже пришли к заключению, что решение спектральных задач — достаточно непростое дело. И это действительно так.Это одна из немногочисленных тем в данной книге, в которой мы несоветуем вам избегать готовых решений».Большинство использованных до настоящего времени пакетов длярешения спектральных задач восходят корнями к монографии [A5].Реализация методов, изложенных в этой книге, получила широкое распространение в виде открытого набора програм EISPACK на Фортране.

Более современной библиотекой, пришедшей ей на смену, сталаLAPACK, также написанная на Фортране с использованием подпрограмм линейной алгебры LINPACK. Библиотека LAPACK отличается более высокой производительностью на современных компьютерах,портирована на Си (CLAPACK), имеет расширения для параллельныхвычислений (ScaLAPACK).Упражнения1) Используя метод прогонки для решения систем линейных уравнений с трёхдиагональными матрицами, напишите программныйкод для интерполяции таблицы значений кубическими сплайнами [1, с.

83]. Исследуйте зависимость погрешности аппроксимациифункции Рунге от количества узлов равномерной сетки, сравните с результатом полиномиальной интерполяции — см. задачу 7на с. 85 [1].2) Найти уровень энергии и волновую функцию 0 () основного состояния в потенциальной яме (), решая конечномерный аналог27спектральной задачи для одномерного стационарного уравненияШрёдингера(︂)︂1 2+ () − (, ) = 0,−|()| → 0 при → 0.2 2ˆ = 0 трёхДля поиска наименьшего собственного значения ˆдиагональной матрицы использовать метод обратных итераций. Проверить работу программы, сравнив с точным решениемдля () = 21 2 .3) В условиях предыдущей задачи найти волновую функцию и уровень энергии первого возбуждённого состояния.4) Что получится, если при поиске первого возбуждённого состояния методом обратных итераций выполнять ортогонализацию с|0 ⟩ (37) только перед началом итерационного процесса? Толькопо окончании итерационного процесса? Коммутируют ли итерации (33) с вычислением проекции (37)? Как это соотносится с темˆ , и проектор на ортогональное дополнефактом, что и матрица ние к |0 ⟩ имеют диагональный вид в базисе |0 ⟩, .

. . , | ⟩?5) Пусть собственные значения матрицы удовлетворяют условию|1 | > |2 | > |3 | > . . ., причем |1 /2 | = |2 /3 | = . . . = 1,1.Считая, что для нахождения каждого из собственных значенийиспользовалось 100 итераций степенного метода (33), оцените погрешность вычисления 1 , 2 , . . ., полагая, что в начальном приближении (32) 1 ≈ 2 ≈ . . .6) Пусть требуется найти наибольшее по модулю собственное значение 1 матрицы порядка = 500 с точностью = 10−4 . Известно, что |1 |/|2 | ≈ 1,01.

Что быстрее: использовать степеннойметод (33) непосредственно для матрицы или предварительно вычислить 2 ? Как изменится ответ, если требуется точность = 10−12 ? Если при этом размер матрицы = 1000?7) Как зависит число арифметических операций на каждом шагеметода вращений Якоби от размера матрицы ? Можно ли сократить число операций в () раз на каждом шаге метода (заисключением первого), используя тот факт, что на каждой итерации меняются только две строки и два столбца матрицы ?8) Получите априорную оценку погрешности вычисления собственных значений в методе Якоби для симметричной матрицы × 28после вращений.