Федченко Г.М. Лекции. Численные методы

Подождите немного. Документ загружается.

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ

Учитывая цикличность выполняемых операций и нецелесообразность

хранения промежуточных результатов, можно записать простой алгоритм

решения СЛАУ методом Гаусса:

1. для k = 1, …, n–1

2. для i = k+1, …, n

5. для j = k+1, …, n

8. для k = n–1, …, 1

Строки 1-6 реализуют прямой ход метода Гаусса, а строки 7-9 – обратный.

Множители стоит сохранять, если система будет решаться несколько раз с

различными правыми частями. Кстати их можно сохранить на

освободившихся местах матрицы A:

Подсчитаем количество операций умножения и деления (как наиболее

трудоемких), выполняемых в ходе работы алгоритма.

Прямой ход. При вычислении выполняется 1+2+…+(n–1)= делений, столько

же умножений при вычислении . пересчитываются (n–1)+(n–2)+…+1 =



раз.

Обратный ход. При вычислении выполняется 1+2+…+(n–1) умножений и n

делений: операций.

Всего операций типа умножения.

Так как при выполнении машинных вычислений неизбежны ошибки

округления, то могут возникнуть ситуации, когда на месте ведущего

элемента окажется 0 из-за потери значимости, что приведет к ошибке

«деление на ноль», или очень маленькое по модулю число, что приведет к

катастрофическому нарастанию погрешности и неверным результатам.

Для того, чтобы избавиться от указанных недостатков классического

метода Гаусса и повысить устойчивость алгоритма к погрешностям

исходных данных и результатов промежуточных вычислений используют

схемы Гаусса с выбором главного элемента на каждом этапе прямого хода.

Существует три стратегии вы бора главного элемента:

1. Выбор главного элемента по столбцу: на каждом этапе прямого хода

сначала в первом столбце промежуточной матрицы, т.е. среди

элементов (i=k, …,n) выбирают наибольший по модулю (главный)

элемент и меняют местами две строки матрицы с тем, чтобы этот



+…+n

=n(n+1)(2n+1)/6

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ

элемент стал ведущим. Таким образом деление будет производиться на

наибольший в данном столбце элемент.

Частичное упорядочивание по столбцам требует внесения в алгоритм

метода Гаусса изменений между строками 1 и 2:

1.1. найти mk такое, что

1.2. если = 0

1.3. то стоп: однозначного решения системы нет

1.4. иначе

1.5.

1.6. для j = k, …, n

1.7.

Разумнее сравнивать не с нулем, а с некоторым малым >0.

Вычисления останавливаются или берутся под особый контроль, если

<. Заметим также, что строка 1.2 позволяет в процессе выполнения

алгоритма проводить исследование системы (5.1) на однозначную

разрешимость.

2. Выбор главного элемента по строке: на каждом этапе прямого хода

ищется максимальный по модулю элемент в строке ведущего элемента

и переставляются местами столбцы. Но перестановка столбцов

означает переобозначение неизвестных, поэтому требуется вводить

дополнительный массив индексов, подверженный аналогичным

перестановкам с тем, чтобы выполнить обратную замену.

3. Выбор главного элемента по всей матрице: на каждом этапе прямого

хода ищется максимальный по модулю элемент среди всех элементов

промежуточной матрицы. Перестановки также требуется запоминать.

Вторая и третья стратегии безусловно сильно усложняют алгоритм.

5.4. Применение метода Гаусса к вычислению определителя и

обращению матрицы

В результате выполнения прямого ходя классического метода Гаусса

матрица A приводится к треугольному виду преобразованиями, не

изменяющими ее определителя. А определитель треугольной матрицы равен

произведению ее диагональных элементов:

detA  = . (5.7)

Таким образом, определитель матрицы равен произведению всех

ведущих элементов при ее преобразовании методом Гаусса.

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ

Для получения определителя матрицы A дополнительно к решению

СЛАУ алгоритм предыдущего пункта должен быть дополнен одной строкой:

10.

Замечание 1. Если алгоритм использовать только для вычисления

определителя, то необходимо изъять из алгоритма строки 4, 7-9.

Замечание 2. Так как перестановка строк или столбцов матрицы меняет знак

определителя, то при использовании метода Гаусса с выбором главного

элемента надо учесть число p произведенных перестановок, а точнее,

четность этого числа, Таким образом, при использовании алгоритма Гаусса с

частичным упорядочиванием строка 10 будет выглядеть следующим

образом:

10.

Для получения матрицы A, обратной к матрице A=, будем исходить

из того, что она является решением матричного уравнения

AX=E. (5.8)

Представляя искомую матрицу X=как набор векторов-столбцов

, , …, ,

а единичную матрицу E как набор единичных векторов

, , …, ,

матричное уравнение (5.8) заменим эквивалентной системой векторно-

матричных уравнений

, , …, . (5.9)

Каждое из уравнений может быть решено методом Гаусса, причем все

они имеют одну и ту же матрицу коэффициентов. Следовательно часть

алгоритма – приведение матрицы A к треугольному виду – общая для всех

систем. Алгоритм можно скорректировать для вычисления обратной

матрицы, изменив строки 4, 7, 9 так, чтобы в роли столбца свободных членов

b оказались единичные векторы:

4.1. для j = 1, …, n

4.2.

7.1. для j = 1, …, n

7.2.

8. для k = n–1, …, 1

9.

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ

В ходе работы алгоритма столбец за столбцом, снизу вверх будет

заполняться обратная матрица X=A.

5.5. Решение линейных систем с помощью LU-разложения

Метод основывается на представлении матрицы A в виде

произведения нижней и верхней треугольных матриц.

Теорема 5.1. Если все главные миноры квадратной матрицы A отличны от

нуля, то существуют такие нижняя L и верхняя U треугольные матрицы, что

A = LU. Если элементы диагонали одной из матриц фиксированы

(ненулевые), то такое разложение единственно.

Напомним, что главными минорами матрицы A=называются

определители подматриц при k=1, 2, …, n–1.

Теорема 5.2. Произвольная невырожденная матрица A перестановкой строк

(столбцов) может быть приведена к матрице с отличными от нуля главными

минорами.

Зафиксируем диагональные элементы матрицы L равными 1. Тогда

необходимо найти такие элементы (при i>j) и (при ij) такие, что

Перемножим матрицы и поэлементно приравняем левые и правые

части:

…

………… ………………… … ……………………

…

Получили nn уравнений относительно nn переменных

Все эти элементы матриц L и U, отличные от 0 и 1, удобно хранить в виде

такой матрицы, равной L–E+U.

Эти неизвестные можно найти в определенном порядке, указанном

стрелками: слева направо и сверху вниз.

Из первой строки уравнений имеем

(j=1, …,n);

из оставшейся части первого столбца уравнений:

(i=2, …,n);

из оставшейся части второй строки:

(j=2, …,n);

из оставшейся части второго столбца:

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ

(i=3, …,n)

и т.д. Последним находится элемент

Общий вид формул для вычисления отличных от 0 и 1 элементов

матриц L и U:

, ij,

(5.10)

, ij.

Организация вычислений по этим формулам должна предусматривать

переключение счета с одной формулы на другую в соответствии с процессом

получения неизвестных, приведших к этим формулам.

Получив LU-разложение матрицы A, можно записать эквивалентное (5.1)

уравнение

LUx=b, (5.11)

решение которого сводится к решению системы

. (5.12)

Таким образом, решение исходной системы сводится к последовательному

решению двух систем с треугольными матрицами коэффициентов.

1. Ly = b:

Следовательно, y при i=1, 2, …, n могут быть найдены по формулам:

. (5.13)

2. Ux = y:

откуда и вычисляются значения неизвестных x в обратном порядке i=n, n–1,

…, 1 по формулам:

. (5.14)

Итак, решение СЛАУ посредством LU-разложения матрицы A

сводится к организации вычислений по формулам (5.10), (5.13), (5.14). На

самом деле результат вычислений по формулам (5.10), (5.13) соответствует

результату прямого хода классического метода Гаусса, а (5.14) – обратного.

Таким образом решение системы LU-разложением – это просто другая схема

вычислений метода Гаусса. Она называется компактной схемой Гаусса или

схемой Холецкого.

5.6. Контрольные вопросы и упражнения

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ

1. Оцените время, затрачиваемое современным персональным

компьютером для решения системы 100, 1000 уравнений методами

Крамера, Гаусса.

2. Подсчитайте количество операций типа умножения, необходимых для

вычисления определителя методом Гаусса.

3. Подсчитайте количество операций типа умножения, необходимых для

вычисления обратной матрицы методом Гаусса.

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ (ПРОДОЛЖЕНИЕ)

Лекция 6

РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ.

ПРЯМЫЕ МЕТОДЫ (ПРОДОЛЖЕНИЕ)

6.1. Метод квадратных корней

Объем вычислений, требующихся для решения линейных

алгебраических задач с симметричными матрицами, можно сократить почти

вдвое, если учитывать симметрию при треугольном разложении.

Пусть A= – симметричная матрица, т.е. . Построим ее разложение в

виде A=UU, где

, .

Перемножим матрицы, получим n(n+1)/2 уравнений относительно

такого же количества переменных (элементов матрицы U):

…

… ……………………

Из первой строки уравнений находим

, (j=2, …,n).

Из второй строки:

, (j=3, …,n).

И наконец:

Общий вид формул для вычисления ненулевых элементов матриц U и U:

, i=1, …, n;

(6.1)

, j=2, …,n.

Получив UU-разложение симметричной матрицы A, решение уравнения

UUx=b

сводится к решению системы

. (6.2)

Таким образом, решение исходной системы сводится к последовательному

решению двух систем с треугольными матрицами коэффициентов.

1. Uy = b:

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ (ПРОДОЛЖЕНИЕ)

Следовательно, y при i=1, 2, …, n могут быть найдены по формулам:

. (6.3)

2. Ux = y:

откуда и вычисляются значения неизвестных x в обратном порядке i=n, n–1,

…, 1 по формулам:

. (6.4)

Решение СЛАУ посредством UU-разложения матрицы A по формулам

(6.3)–(6.5) называется методом квадратных корней или схемой Холецкого.

6.2. Метод прогонки решения систем с трехдиагональными

матрицами коэффициентов

Часто возникает необходимость в решении линейных алгебраических

систем, матрицы которых являются слабо заполненными, т.е. содержащими

много нулевых элементов, а ненулевые элементы выстраиваются в

определенную структуру. Среди таких систем выделим системы с матрицами

ленточной структуры, в которых ненулевые элементы располагаются на

главной диагонали и на нескольких побочных диагоналях. Для решения

систем с ленточными матрицами коэффициентов метод Гаусса можно

трансформировать в более эффективные методы.

Рассмотрим случай системы линейных уравнений с матрицей

коэффициентов ленточной структуры шириной три. К таким системам

приводят задачи о сплайн-интерполяции, о решении разностными методами

обыкновенных дифференциальных уравнений и уравнений в частных

производных. Система линейных уравнений имеет трехдиагональную

структуру:

, (6.5)

т.е. каждое уравнение связывает три соседних неизвестных:

, i=1, 2, …, n,

. (6.6)

Выразим из первого уравнения x через x:

;

из второго x через x:

 

;

и т.д. Из последнего уравнения:

 

Таким образом, для i=1, 2, …, n:

, где (6.7)

, . (6.8)

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ (ПРОДОЛЖЕНИЕ)

Так как , и .

Таким образом, решение систему линейных уравнений (6.6) сводится

к двум этапам:

1. Вычисление прогоночных коэффициентов P, Q, i=1, 2, …, n по

формулам (6.8) – прямая прогонка.

2. Вычисление неизвестных x при i=n, n–1, …, 1 по формулам (6.7)

– обратная прогонка.

Замечание 1. Обозначим , тогда , при i=2,…,n.

Непосредственной проверкой можно убедиться, что имеет место

представление A=LU, где

Таким образом, LU-разложение трехдиагональной матрицы A может быть

выполнено простым алгоритмом, вычисляющим R, P при возрастающих

значениях i.

Замечание 2. Можно вычислить и определитель матрицы A

det A= . (6.9)

6.3. Контроль точности и уточнение приближенного решения в

рамках прямых методов

Прямые методы приводят к точному решению СЛАУ при точном

выполнении предусматриваемых соответствующими алгоритмами

арифметических операций. Реальные же вычисления производятся

приближенно. Как отражается на результате решения системы подмена

арифметики действительных чисел машинной арифметикой, зависит от

самой решаемой системы, параметров применяемого компьютера, способов

реализации алгоритма. В любом случае, практически вместо точного

решения СЛАУ прямой метод дает приближенное решение.

Обозначим полученное приближенное решение через и подставим в

исходное уравнение. Вектор называют невязкой приближенного решения .

По величине невязки можно с некоторой осторожностью судить о близости

найденного решения к точному решению x. Если |||| – недостаточно малая

величина, то следует искать вектор поправку такой, что , т.е.

Следовательно

Нахождение поправки сводится к решению той же системы

уравнений, где в качестве столбца свободных членов берется вектор невязки.

РЕШЕНИЕ СЛАУ. ПРЯМЫЕ МЕТОДЫ (ПРОДОЛЖЕНИЕ)

Так как матрица коэффициентов левой части системы не меняется, нет

надобности прямого хода преобразований коэффициентов при неизвестных

(другими словами – LU-разложения), достаточно выполнить действия с

новыми свободными членами.

Прибавив найденную поправку, получаем уточненное решение . Если

величина |||| (или ||||/ ||||) окажется недостаточно малой, процесс уточнения

может быть повторен: ищется поправка как приближенное решение системы

, где ; тогда более точным должно быть решение .

Но сходимость к нулю невязок в таком процессе может и не

наблюдаться. Обычно делают не более двух-трех шагов уточнения, причем

рекомендуется производить вычисление невязок в режиме накопления. если

в этом процессе не происходит сближения при k=2,3,…, то это говорит

скорее всего о плохой обусловленности данной системы и о том, что ее

решение не может быть уточнено таким способом.

Описанный контроль точности по невязкам не требует большого

увеличения вычислительных затрат, но требуется отводить дополнительное

место в памяти под хранение исходных данных.

6.4. Обусловленность матриц

Не всегда малость невязки гарантирует малость погрешности

найденного решения.

Замечание 3. Гауссово исключение с выбором главного элемента

гарантированно дает малые невязки. Малость здесь трактуется относительно,

т.е. по отношению к элементам матрицы A, элементам матриц в

промежуточных вычислениях и к компонентам вычисленного решения.



||||x||∙||A||∙



Но даже если невязка мала, вектор ошибки не обязательно мал:

||x||||x||∙cond(A)∙



Обусловленность – это внутреннее свойство матрицы A, не зависит от

метода решения СЛАУ. Матрицы с большим числом обусловленности дают

большие ошибки при решении систем. Логарифм числа обусловленности

приближенно равен числу значащих цифр, теряемых в решении системы

Ax=b.

Например, при



10

–7

, cond(A)10

, относительная ошибка решения будет

10

–3

, т.е. рассчитывать можно на три верных разряда в решении.

Число cond(A) измеряет насколько близка к вырожденной матрица A

и, что еще важнее, насколько чувствительно решение системы Ax=b к

изменениям в A и b.