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

Как следует из определения,cond() ≥ 1 для любой невырожденной матрицы ; в случае, если матрица вырожденная, знаменатель (15) обращается в ноль.3 Здесь предполагается, что произведения −1 не настолько малы, чтобы происходила потеря порядка.12Для эрмитовых матриц выражение (15) может быть записано в видеcond() = max | |/ min | |, где — собственные числа.В полном соответствии с рассмотренным выше примером, из формул (14) и (15) видно, что погрешность численного решения системы оказывается нечувствительной к масштабным преобразованиям ← −1 .

Также из (15) следует, что повышенная погрешностьчисленного решения системы связана с наличием в матрице одновременно больших и малых масштабов. Сложение чисел различных порядков величины приводит к отбрасыванию части значащих разрядовменьшего по модулю слагаемого (вплоть до катастрофической потериточности, рассмотренной в п.

1.2). Поскольку количество разрядов, отбрасываемых при сложении + , приблизительно равно модулю разности порядков и , т. е. | log |/||, можно ожидать, что логарифмчисла обусловленности (15) будет приближённо равен числу теряемыхзначащих цифр при решении системы (1).1.4. Определитель и обратная матрицаРассмотренный выше метод исключения Гаусса не только позволяет решать системы линейных уравнений, но и является основой дляэффективных способов вычисления определителя и обратной матрицы.

Обратим внимание, что преобразование (3) квадратной матрицы оставляет неизменным её определитель. Полученная в результатепрямого хода метода исключения матрица имеет верхнетреугольныйвид (4), и её определитель равен произведению диагональных элементов. (В этом несложно убедиться, вычисляя определитель матрицы (4)разложением по первому столбцу.) Число операций, необходимых длявычисления определителя рассмотренным способом, приблизительноравно 23 3 + ≈ 23 3 . При вычислении определителей матриц большогоразмера и в случае, если матричные коэффициенты сильно отличаются по модулю от 1, следует иметь в виду возможность переполнения ипотери порядка.Рассмотрим задачу о нахождении обратной матрицы −1 , полагая, как и раньше, невырожденной квадратной матрицей × .

Аналогично формуле Крамера (2) для решения линейных систем и нахождению определителя матрицы через разложение по строке (столбцу),хорошо известное из алгебры выражение для обратной матрицы−1 =adj()det являет собой ещё один пример неэффективного способа вычисления13(здесь adj() — присоединённая матрица, составленная из алгебраических дополнений). Значительно более быстрый и притом относительнонесложный в реализации способ состоит в следующем.

По определениюобратной матрицы (−1 ) = ,(16)где — символ Кронекера; выражение (16) при фиксированном можно рассматривать как систему линейных алгебраических уравненийна неизвестных — матричные элементы -го столбца матрицы −1 .Несмотря на то, что формально для нахождения всех матричных коэффициентов −1 требуется решить линейных систем (16) — по одной для каждого значения = 1, .

. . , , — необходимое для этого числоарифметических операций при правильной организации вычисленийвозрастёт лишь приблизительно в 2,5 раза, оставаясь равным (3 ).Действительно, заметим, что преобразование (3) матрицы во времяпрямого хода метода исключения Гаусса никак не зависит от вектора b правых частей системы. Поскольку системы уравнений (16) при = 1, . . . , имеют общую квадратную матрицу , отличаясь лишьвектором правых частей уравнений b1 , . .

. , bn , можно обобщить выражение (3) для эффективного решения систем (16): ← / ;← 0;← − · , = + 1, . . . , ;← − · , = + 1, . . . , ,(17)где — квадратная матрица, составленная из векторов в правой части систем (16) при = 1, . . . , ; её начальное значение = .1.5. Метод прогонки для трёхдиагональных матрицВ предыдущих пунктах мы рассматривали общий случай системылинейных уравнений с неизвестными. Вместе с тем в физике есть множество частных случаев, в которых решение систем уравнений специального вида может быть выполнено значительно эффективнее.

С одним из примеров таких систем мы уже знакомились при рассмотрении задачи полиномиальной интерполяции [1, с. 72]. Для определениякоэффициентов полинома () = 0 + 1 1 + . . . + с помощьюрешения системы алгебраических уравнений требуется (3 ) арифметических операций, тогда как при использовании интерполяционныхполиномов Лагранжа или Ньютона аналогичный «подготовительный»14этап вычисления коэффициентов может быть выполнен всего за (2 )операций.Другим чрезвычайно важным примером являются системы уравнений сматрицами, т.

е. матрицами, большинство элементов которых равны нулю. В этом пункте мы познакомимся с хрестоматийным методом прогонки для трёхдиагональных матриц — по сути,с частным случаем рассмотренного выше метода исключения Гаусса.Трёхдиагональные матрицы уже встречались нам при обсуждении интерполяционных кубических сплайнов [1, с. 83]. Кроме того, ленточныематрицы регулярно возникают при решении дифференциальных уравнений в частных производных с использованием конечно-разностных иконечно-элементных методов, знакомство с которыми нам ещё предстоит. Размер таких матриц обычно равен числу узлов сетки (а в некоторых случаях может и превосходить его в несколько раз) и, следовательно, может быть достаточно большим — вплоть до десятков и сотен тысяч.

Использование специальной (разреженной) структуры таких матриц позволяет повысить эффективность расчётов на много порядковвеличины.Запишем систему уравнений с трёхдиагональной матрицей в видеразреженными −1 + + +1 = , = 1, . . . , ,1 = = 0.(18)Прямой ход метода исключения Гаусса (3) для системы (18) сводитсяк исключению поддиагональных элементов , = 2, . . . , :←, ← 0, ← − −1 , ← − −1 .(19)−1Сумма в формуле для обратного хода (5) в случае трёхдиагональнойматрицы будет содержать всего один член:)︀1 (︀; = − +1 , = − 1, . .

. , 1.(20) =Решение системы методом прогонки (19), (20) требует всего 4 ячеек памяти для хранения матричных коэффициентов (18) и 8 арифметических операций, что значительно экономичнее метода исключенияГаусса (3), (5) для матриц общего вида.Можно показать [2, гл. 3], что достаточным условием устойчивостиметода прогонки может служить: | | ≥ | | + | | со строгим неравенством хотя бы для одного .Подчеркнём, что данное условие является достаточным, но не необходимым, и в физических расчётах прогонка в большинстве случаевоказывается достаточно устойчивой даже при нарушении преобладания диагональных элементов.15товпреобладание диагональных элемен-1.6. Модификация метода прогонки для периодических граничных условийРассмотренный в предыдущем пункте метод прогонки будет применён нами в п. 3.3 для решения уравнения теплопроводности − = (, ) с использованием неявных схем.

Заменяя частную производную на конечную разность ℎ12 (+1 − 2 + −1 ), мы получим систему уравнений, каждое из которых связывает значения температурыв трёх соседних узлах. Использование нулевых условий на границе(0, ) = (, ) = 0 обеспечит зануление двух матричных коэффициентов в первом и последнем уравнении, так что мы получим системус трёхдиагональной матрицей (18). Что делать, если нам необходиморешить уравнение теплопроводности на цилиндре, используя периодические граничные условия (0, ) = (2, ), (0, ) = (2, )? Очевидно, что при этом каждое уравнение, в том числе самое первое ипоследнее, содержит по три члена, ни один из которых не равен нулю, так что матрица системы уже не является трёхдиагональной (нидаже ленточной).

Рассмотрим кратко один из экономичных способоврешения таких систем, предложенный в работе [A1]. Запишем исходную систему линейных уравнений в виде(21)x = d,где вектор x — искомое решение, d — известный вектор — матрица системы (21) порядка :⎛1 1 0 . . .001⎜ 2 2 2 . . .000⎜⎜ 0 3 3 . . .000⎜=⎜ ............⎜ ......⎜⎝ 0 0 0 .

. . −1 −1 −1 0 0 . . .0правых частей,⎞⎟⎟⎟⎟⎟⎟⎟⎠(22).Обратим внимание, что первые − 1 строк и столбцов матрицы образуют трёхдиагональную подматрицу — обозначим её ′ . Перепишемисходную систему (21) через ′ , что позволит нам решать возникающиеуравнения методом прогонки:(︂ ′)︂ (︂ ′ )︂ (︂ ′ )︂u′xd=,(23)v′T где обозначает транспонирование, штрихованными буквами обозначены векторы из (−1)-мерного подпространства. В частности, u′ , v′ —16векторы, состоящие из первых ( − 1) элементов -го столбца и -йстроки (22):u′=(10...0 −1 )T ,v′=(0...0 )T ,′d=(12...(24)T−1 ) .Выражая x′ из первого (векторного) уравнения системы (23):−1x′ = (′ )(d′ − u′ ) ≡ p′ + q′ .(25)Здесь ( − 1)-мерные векторы p′ и q′ удовлетворяют системам уравнений с трёхдиагональной матрицей ′ :′ p′ = d′ ,′ q′ = −u′ .(26)Выразим из второго (скалярного) уравнения системы (23): =)︀1 (︀ − (v′ , x′ ) .(27)Подставляя x′ из (25) в (27) и раскрывая скалярное произведение, получаем окончательно для : = − ′1 − ′−1 − (v′ , p′ ).=′ + (v′ , q′ ) + 1′ + −1(28)Таким образом, для решения задачи (21) необходимо выполнитьследующие шаги:1) найти p′ , q′ (26) методом прогонки (19), (20);2) вычислить по формуле (28);3) вычислить x′ = p′ + q′ (25).Вычисления по указанному алгоритму потребуют порядка 14 операций, что почти вдвое больше, чем для «обычного» метода прогонки(19), (20) для трёхдиагональных матриц.

Действительно, на шаге 1дважды используется метод прогонки, однако ввиду того, что обе системы (26) имеют одну и ту же матрицу ′ , их одновременное решениевозможно за 12 операций. Ещё 2 операций требуется на шаге 3.171.7. Спектральная задачаРасчёт резонаторов и волноводов в оптике, определение резонансных частот для классических механических систем, поиск уровнейэнергии квантовых частиц и ряд других численных задач сводятся кпоиску собственных значений и собственных векторов матриц:u = u,(29)где и u — собственное значение и собственный вектор матрицы порядка . Из курса алгебры известно, что необходимым и достаточнымусловием существования нетривиального решения u (29) служит равенство нулю определителя0 = det | − | ≡ (),(30)где — единичная матрица.