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

Быстрый (экспоненциальный) рост числа операций при увеличении размераматрицы делает метод Крамера непригодным для вычислений определителей матриц даже небольшого размера. Как мы вскоре увидим,1 Нелинейные системы могут быть решены с использованием итераций, при этомна 2каждом шагенеобходимо решать систему линейных уравнений.1,25 × 1017 операций в секунду согласно рейтингу www.top500.org.6эффективная стратегия состоит в использовании метода решения системы линейных уравнений для вычисления определителя, а не наоборот, как предполагает формула Крамера (2).Для начала рассмотрим наиболее простую модификацию метода(как будет показано в следующем параграфе, она может давать совершенно неудовлетворительные результаты даже для хорошо обусловленных систем). В данном рассмотрении мы будем использовать длякоэффициентов системы, изменённых в процессе решения, те же обозначения , , что и для их первоначальных (входных) значений.Жертвуя при этом математической строгостью обозначений, взаменмы получаем более тесную связь формул решения с компьютернойпрограммой, в которой содержимое массивов (матриц) изменяется впроцессе решения, в то время как имена массивов остаются прежними.Решение системы уравнений методом исключения Гаусса происходит в два этапа, которые называютиметода.заключается в приведении матрицы системы к верхнетреугольному виду (так, чтобы все элементы ниже главной диагонали были равны нулю: = 0 ∀ > ).

Зануление производится последовательно по столбцам: = 1, 2, . . . , − 1, а в пределах каждого столбца —по строкам: = + 1, + 2, . . . , . Для зануления коэффициента нужно вычислить отношение = / , после чего вычесть из -гоуравнения системы -е, умноженное на :прямым обратным ходомПрямой ход ← / ;← 0;← − · , = + 1, . . .

, ;(3)← − · .Здесь и далее символ ← обозначает оператор присваивания, т. е.запись ← подразумевает, что значение переменной записываетсяв переменную (ячейку памяти) . Повторение процесса (3) в цикле длястрок = +1, . . . , приведёт к занулению всех элементов -го столбцаматрицы ниже главной диагонали. Выполняя описанные выше операции последовательно для столбцов = 1, 2, . . . , − 1, мы преобразуемматрицу системы к верхнетреугольному виду:7⎛⎜⎜⎜⎜⎜⎝110...1222...0000......1,−12,−1....

. . −1,−1...012...−1,⎞⎛12...⎟⎜⎟⎜⎟⎜⎟⎜⎟⎜⎠ ⎝ −1⎞⎛12...⎟ ⎜⎟ ⎜⎟ ⎜⎟=⎜⎟ ⎜⎠ ⎝ −1⎞⎟⎟⎟⎟ . (4)⎟⎠После этого решение системы x может быть легко найдено с помощьюметода исключения Гаусса: вначале вычисляетсяпоследняя компонента = / , затем найденное значение подставляется в предпоследнее уравнение: −1 = (−1 − −1, )/−1,−1и т. д. Общий вид формулы обратного хода:⎞⎛∑︁1 ⎝ ⎠ , = − 1, − 2, . . . , 1.(5) = −=+1обратного ходаПодсчитаем число операций, необходимых для решения системыметодом исключения Гаусса. Вычитание -го уравнения из -го в соответствии с (3) требует ( − + 1) операций умножения, столько жевычитаний и одно деление — для простоты будем считать их вместе как2( − + 1) операцию. Общее число операций прямого хода при этоместь1 =−1∑︁∑︁2( − + 1) = 2=1 =+1−1∑︁( − )( − + 1) ≈=12 3 .3(6)(Данная асимптотика может быть получена переходом к суммированию по ′ = − и последующей заменой дискретной суммой на интеграл.) Число операций обратного хода (5) есть2 =∑︁(1 + 2( − )) = (2 ).(7)=1Поскольку 1 = (3 ), числом операций обратного хода можно пренебречь, так что полное число операций Σ = 1 + 2 ≈ 1 ≈ 32 3 .81.2.

Метод исключения Гаусса с выбором главногоэлементаНапишем программу в соответствии с полученными выше формулами (3), (5) и попробуем решить следующую систему уравнений:⎞⎛⎞ ⎛⎞⎛1−1050,2 1⎝ 1 0,04 1 ⎠ ⎝ 2 ⎠ = ⎝ −10 ⎠ .(8)320−511Полученные в результате корни 1 = 2 = 0, 3 = −10 при подстановке в систему (8) не дадут равенства правой и левой частей. Разберёмся, почему так происходит.На первом шаге прямого хода метода исключения Гаусса первоеуравнение (8) должно быть умножено на = 21 /11 = 1/5 = 0,2 ивычтено из второго уравнения.

Выполняя эти вычисления в уме, мыполучим во второй строке матрицы (0; 0; 0,8). Обратим внимание назануление диагонального элемента 22 = 0: на следующем шаге метода исключения (3) это должно привести к делению на 0 с аварийнымостановом программы или появлению в расчётах нечисловых значенийinf и nan.Однако, поскольку входящие в матрицу (8) десятичные числа 0,2и 0,04 непредставимы в виде конечных двоичных дробей, при вычислении на ЭВМ возникает погрешность округления, в результате чеговместо 22 = 0 мы получим хотя и очень малое, но отличное от нулязначение 22 = −6,94 · 10−18 :⃒⎛⎞50,21 ⃒⃒ −10Λ(1) = ⎝ 0 −6,94 · 10−18 0,8 ⃒⃒ −8 ⎠ .(9)01,22 ⃒ 10(Здесь Λ(1) — расширенная матрица 3 × 4 системы уравнений после зануления столбца = 1.) Поскольку 22 ̸= 0, в дальнейших вычислениях не возникнет деления на ноль — мы не получим ни аварийногоостанова программы, ни значений inf и/или nan, которые могли бысигнализировать о возникновении нештатной ситуации в расчётах.Для зануления элемента 32 в столбце = 2 в соответствии с (3)необходимо вычесть из третьей строки (9) вторую, умноженную на =1,2/(−6,94 · 10−18 ) = −1,73 · 10+17 .

Чтобы понять, чему будет равенрезультат машинной операции32 ← 2 − 0,8 · (−1,73 · 10+17 ) = 2 + 1,38 · 10+17 ,9(10)вспомним, что компьютер оперирует с числами с плавающей точкой(запятой), имея в распоряжении 52 двоичных разряда [1, п. 2.2]. Единица в младшем разряде (ULP, или машинное эпсилон) равно =0,00 . . .

012 = 2−52 ≈ 2,22 · 10−16 и характеризуетточность вычислений. Поскольку операнды в (10) отличаются более чем в раз, можно ожидать, что 32 ← 2+1,38·10+17 = 1,38·10+17 . Подчеркнём, что результат машинного сложения в (10)(в машинном представлении) равен второму слагаемому!Произошедшее является примером: несмотря на то, что результат операции (10) имеет все 52 верных двоичных разряда, при вычитании строки 2 матрицы (9) из строки3 (3 ← 3 − 2 ) информация о третьей строке была полностьюутрачена, оказавшись за пределами точности вычислений:⃒⎛⎞⃒50,21−10⃒⃒⎠.0,8−8Λ(2) = ⎝ 0 −6,94 · 10−18(11)⃒001,38 · 10+17 ⃒ −1,38 · 10+18относительнуюв точности побитовокатастрофической потери точ-ностиКак следствие, строки 2 и 3 матрицы (11) оказались почти пропорциональны друг другу.

Используя формулы обратного хода (5), несложнополучить из (11) ответ x = (0; 0; −10), который, очевидно, не удовлетворяет третьему уравнению исходной системы (8).Заметим, что, хотя коэффициенты в (9) и даже в (8) уже содержатпогрешность округления, это оказывает лишь незначительное влияниена ответ: прямым вычислением можно убедиться, что относительнаяпогрешность вычисления корней (8) и (9) имеет тот же порядок величины, что и погрешность округления . Причиной катастрофическойпотери точности является вычитание двух строк матрицы с домножением на большой коэффициент (3), что, в свою очередь, связано свозникновением в процессе решения близкого к нулю диагональногоэлемента матрицы системы.Чтобы избежать рассмотренной выше потери точности, достаточно лишь немного модифицировать метод исключения Гаусса, дополнивегона каждом шаге прямого хода.

Аименно, перед применением формулы (3) при фиксированном и всех = +1, . . . , необходимо вначале найти строку с максимальным помодулю элементом и поменять её местами с -й строкой расширенной матрицы системы. В программной реализации метода более простой и эффективной альтернативой перестановке строк с копированиемданных может служить массив перестановок индексов, через которыйбудет вестись обращение к массиву матричных элементов: например,выбором главного элемента10вместо A[i][j]-=c*A[k][j] формула (3) может быть запрограммирована как A[p[i]][j]-=c*A[p[k]][j], где A — расширенная матрица системы, int p[n] — массив перестановок индексов. Массив p[n] вначаледолжен быть инициализирован единичной перестановкой p[i]=i, а перед каждым шагом прямого хода необходимо менять местами значения в ячейках p[k] ↔ p[m]. Видно, что выбор главного элемента лишьнемного усложняет код программы, но существенно повышает её надёжность, позволяя предотвратить потерю точности.В заключение данного параграфа ещё раз обратим внимание на то,как важно помнить о специфике машинных вычислений и критическиотноситься к результатам, выдаваемым даже написанной «без ошибок»компьютерной программой.1.3.

Погрешность и невязкаГоворя о любом численном методе, всегда следует иметь в виду двавопроса: его эффективность (количество операций для получения ответа) и точность, которую он обеспечивает. Первый аспект уже былрассмотрен нами выше, теперь поговорим о погрешности численногорешения систем линейных уравнений.Обычно используют две меры погрешности. Первая, наиболее естественная и представляющая наибольший интерес, —,определяемый как разность двух решений: точного (x) и полученного в численных расчётах (x* ):вектор ошибки = x − x* .невязка(12)Вторая —, полученная при подстановке численного решения висходную систему уравнений:r = b − x* = (x − x* ) = · .(13)Вектор невязки (13) может быть легко вычислен непосредственно послеполучения решения системы.

Важным достоинством метода исключения Гаусса с выбором главного элемента является гарантированно малая величина невязки: |r| ≈ || × |x| × , т. е. компоненты вектораневязки (13) по порядку величины будут равны произведению матричных элементов , компонент вектора решений x и машинного .Однако малость вектора невязки r (13) ещё не означает малостипогрешности (12) численного решения. Действительно, как следуетиз (13), = −1 r, т. е. в случае, если матрица является «почтивырожденной», численное решение может иметь большую погрешность|| даже при относительно малой величине невязки |r|.11Какую матрицу следует считать «почти вырожденной»? Первое,что приходит в голову — это малость определителя det . Однако обратим внимание, что если выполнить масштабное преобразование (1),разделив каждый элемент матрицы A на одну и ту же константу ≫ 1 : ← −1 , определитель матрицы уменьшится в большое число раз, при этом вряд ли можно ожидать сколько-нибудьзаметного изменения корней x* и относительной погрешности их вычисления3 .

В самом деле, ЭВМ оперирует числами с плавающей точкой, сохраняя отдельно порядок и мантиссу каждого числа. Умножение всех чисел в расчётах на одну и ту же константу −1 ≪ 1 приведётлишь к уменьшению порядка всех чисел, что никак не отразится наарифметических операциях над матричными элементами . Такимобразом, понятие «почти вырожденной» матрицы не сводится к малости определителя det .Чтобы лучше понять суть дела, вспомним, что при изучении численных методов интегрирования обыкновенных дифференциальных уравнений [1, с.

89] мы уже встречались с примерамизадач, небольшое изменение условий которых приводит к значительному изменению ответа. Наличие малой ошибки (/2) во входных параметрах и погрешностей округления при выполнении промежуточныхвычислений приведут к погрешности численного решения . При одинаковой относительной погрешности входных условий /2 погрешностьответа может существенно отличаться дляизадач.Аналогичное справедливо и в отношении численного решения систем линейных уравнений. Абсолютную погрешность (12) решенияможно оценить как|| ≈ |x* | · cond() · ,(14)плохо обусловленныххорошо плохо обуслов-ленныхгде cond() —число обусловленности матрицы :(︂cond() =maxx|x||x|)︂ ⧸︁ (︂)︂|x|min.x|x|(15)Другими словами, число обусловленности (15) есть отношение максимального и минимального числа раз, в которое матрица изменяет норму произвольного вектора x.