Численные методы линейной алгебры.

Лекция 3. Стр. 1 Лекция 3 Численные методы линейной алгебры Численные методы линейной алгебры: прямые и итерационные методы решения систем линейных алгебраических уравнений (метод Гаусса, метод простой итерации, метод Зейделя). В классических курсах методов вычислений обычно выделяют четыре основные задачи линейной алгебры: решение системы линейных алгебраических уравнений, вычисление определителей, вычисление обратной матрицы, определение собственных значений и собственных векторов матрицы. Мы рассмотрим только первую и последнюю из этих основных задач. Вторую и третью задачи вкратце обсудим при изучении задачи решения линейных систем. Решение систем линейных алгебраических уравнений. Постановка задачи Рассмотрим систему n линейных алгебраических уравнений относительно n неизвестных x1 ,x2 ,...,xn : a11 x1  a12 x2  ...  a1n xn  b1 , a x  a x  ...  a x  b ,  21 1 22 2 2n n 2  ... a n1 x1  a n 2 x2  ...  a nn xn  bn . Решением системы называется совокупность значений неизвестных x1  x1* , x2  x2* , ..., xn  xn* , при подстановке которых все уравнения системы обращаются в тождества. Система, имеющая хотя бы одно решение, называется совместной; система, не имеющая ни одного решения — несовместной. n Линейная система “свернутом” виде может быть записана в виде a j 1 ij x j  bi , i  1, n , или в матричном виде  a11 a12 ... a1n   x1   b1         a 21 a 22 ... a 2 n   x2   b2  , b   , x   , Ax  b , где A    ... ... ... ... ... ...       a       n1 a n 2 ... a nn   xn   bn  матрица A называется матрицей системы, вектор b называется правой частью системы, вектор x — решение системы. Если det A  0 , то матрица системы невырождена, у нее существует обратная матрица. Известно, что в этом случае решение системы существует, единственно и устойчиво по входным данным. Это означает, что рассматриваемая задача корректна. Пусть x   ( x 1 , x  2 , ..., x  n ) - приближенное решение системы. Мы будем стремиться к получению решения, для которого погрешность e  x  x  мала. (Количественные характеристики «величины» погрешности будут введены позднее.) Качество полученного решения далеко не всегда характеризуется тем, насколько мала погрешность e  x  x  . Иногда вполне удо влетворительным является критерий малости вектора невязки r  b  Ax . Вектор r показывает, насколько отличается правая часть системы от левой, если внее подставить приближенное решение. Заметим, что r  b  Ax   Ax  Ax   A( x  x  ) и поэтому погрешность и невязка  1 связаны равенством e  x  x  A r .Решение системы формально легко получить, умножив обе части уравнения Ax  b слева на матрицу A 1 имеем A 1  (Ax)  A 1  b, A 1  A  E, , Ex  x, x  A 1  b Оба эти случая хорошо иллюстрируются геометрически для n  2 : Лекция 3. Стр. 2 Обусловленность задачи решения линейной системы. Приведенные выше рисунки иллюстрируют некорректную задачу и хорошо обусловленную корректную задачу. Однако, система Ax  b эквивалентна системе Ax  b  0 , а поскольку из-за округления чисел при вводе в память компьютера левая часть будет отлична от точного нуля, то можно утверждать, что, вообще говоря, любая система после ввода в память компьютера становится несовместной. Б. Парлетт в книге "Симметричная проблема собственных значений" пишет, предупреждает: "Причудливые подчас эффекты арифметики конечной точности должны быть всегда в поле зрения при матричных вычислениях". Попытаемся понять, как сильно могут влиять погрешности округления на погрешность конечного результата, какие свойства матрицы системы определяют ее "чувствительность" к погрешностям округления и какие "причудливые эффекты" могут возникать в машинных расчетах, т.е. попытаемся исследовать обусловленность задачи решения линейной системы. Обусловленность решения линейной системы важна еще и потому, что матрица системы и правые части в большинстве практических научных расчетах содержат ошибки, мера влияния которых на погрешность результата также определяется обусловленностью задачи. Ниже приведена геометрическая иллюстрация плохо обусловленной задачи. Прямая  получена как результат возмущения правой части уравнения, описывающего прямую  . Малая погрешность правой части приводит к катастрофической погрешности: точка А — решение исходной системы, точка В — решение возмущенной системы. Рассмотрим пример x 2  2, x 2  2,  x1   x1  и    x1  1.0001 x 2  2.0001  x1  1.0001 x 2  2.0002 . Очевидно, что решением первой системы является вектор 1.0000, 1.0000 , а решением второй — 0.0000, 2.0000 , т.е. изменение правой части в пятой значащей цифре привело к изменению решения в первой значащей цифре. Заметим, что определители матриц малы ( 10 4 ). Рассмотрим еще несколько примеров. 6 6 6  2  10 x1  10 x 2  3  10 , 12 Определитель матрицы системы  равен 10 , хотя ее реше6 6 6   10 x1  10 x 2  2  10 ния мало зависят от погрешностей исходных данных. В этом легко убедиться, умножив оба ее уравнения на 10 6  2 x1  x 2  3,  x1  x 2  2. :  С другой стороны, определители матриц систем Лекция 3. Стр. 3  x1  100 x 2  100 ,  x1  100 x 2  100 , и  равны единице.  x2  1 x2  0   Решение первой системы —( 0, 1 , а решение второй системы — 100, 0 . Т.е малая по- 100 2  1  1.00005 , относительная погрешность составляет 100 100 2  1 1%, приводит к катастрофической погрешности решения  x *   100 .00 , отно1   грешность правой части  b   сительная погрешность решения равна 100%. Таким образом, мы убедились, что det A (определитель матрицы системы) не может служить мерой обусловленности матрицы. Норма матрицы. Число обусловленности матрицы. Если в пространстве векторов x  ( x1 , x2 , ..., xn ) введена норма x и если 0  m  Ax x  M , x  0 , числа m и M — наименьшее и наибольшее значения своеобразного "коэффициента растяжения". Ax Согласованной нормой в пространстве матриц называют величину A  max x0 x M. Например, если x  x 1  n n  xi , то согласованная норма матрицы A  A 1  max  aij , j i 1 n если x  x 2  x если x  x e  x i 1 i 1 i   2 , то согласованная норма матрицы A  A 2  max A  A T , 2 , то согласованная норма матрицы A  A e  i n i 1 n i 1 j 1 если x  x i  max xi , то согласованная норма матрицы A  A i  max i i n  a 2 ij , n a j 1 ij . Внешне столь различные, эти нормы эквивалентны: когда некоторая последовательность векторов по одной из этих норм стремится к нулю, то она стремится к нулю и по остальным нормам: lim x ( n ) n N1  0 тогда и только тогда, когда lim x ( n) n ществуют числа  1  0 и  2  0 такие, что  1 A N1  A N2 N2  0 . Или, что то же самое, су- 2 A N1 . Величина m тоже играет важную роль в матричных вычислениях. Если матрица A невыро1 ждена, то она обратима и A  1 . m Абсолютная и относительная погрешности вектора. Абсолютная погрешность вектора x ( x  )  x  x  Относительная погрешность вектора x  ( x )    x  x x . Выбор той или иной конкретной нормы в практических задачах диктуется тем, какие требования предъявляются к точности решения. Пусть x — точное решение линейной системы Ax  b , а x*  x  x — решение системы с возмущенной матрицей A*  A  A : Ax  b , A  Ax  x  b , Лекция 3. Стр. 4 вычитая из второго уравнения первое, имеем Ax  Ax  x  0 , откуда x  A 1 Ax  x  и, следовательно, Величину condA   A  A 1  x x  x  A  A 1  A A  condA  A A . M принято называть числом обусловленности матрицы. m Оценка числа обусловленности. Замечательно, что число обусловленности матрицы совершенно также определяет влияние на погрешности решения возмущения правой части системы. Пусть x*  x  x — решение системы с возмущенной правой частью: b*  b  b : Ax  b , Ax  x  b  b , Ax  b , x  A 1 b , m x  b  Ax  M x и x x   A 1 b x  x b b b M b  condA   A  A 1  condA  т.е. x b m b b b Заметим, что если x  x 2 , то condA   A  A 1       max A  A T . Такое определение  min A  A T числа обусловленности позволяет дать следующую геометрическую интерпретацию обусловленности матрицы. Матрица A отображает n-мерную единичную сферу (множество векторов единичной длины x : x  1 ) в множество векторов различной длины Ax : x  1 . Это мно-     жество — r-мерный эллипсоид в R , r — ранг матрицы А, для невырожденных систем r = n. На рисунке изображен случай n =2. n     Здесь  max A  A T и  min A  A T — наибольшее и наименьшее сингулярные числа матрицы — соответственно большая и малая полуоси эллипса и плохая обусловленность отвечает сильно вытянутым эллипсам. Таким образом, чем больше число обусловленности, тем чувствительнее решение линейной системы к погрешностям округления. Системы (матрицы) с большим числом обусловленности называют плохо обусловленными. Реальное вычисление числа обусловленности предполагает знание нормы обратной матрицы, вычисление которой сравнимо по количеству операций с решением системы. Однако для анализа обусловленности системы достаточно знать какую-либо хорошую оценку числа обусловленности, например, часто вычисляется очень точная оценка, нижняя граница истинного  j  числа обусловленности condA   max A j z y , где A  j  — j-й столбец матрицы А, A T y  e , Az  y , e — специально выбранный вектор, координаты которого  1 . Решение линейной системы методом Гаусса Методы решения систем линейных алгебраических уравнений можно разделить на точные и приближенные. Примером точных методов может служить метод Крамера, который на самом Лекция 3. Стр. 5 деле никогда не применяется из-за огромного количества вычислений и связанных с ними ошибок округления. Точные методы решения линейных систем применяют для решения линейных систем относительно небольшой размерности (до 10 3 ). Метод Гаусса — точный метод решения невырожденной системы линейных алгебраических уравнений. Метод Гаусса, его еще называют методом гауссовых исключений, состоит в том, что систему n линейных алгебраических уравнений относительно n неизвестных x1 , x2 , ..., xn a11 x1  a12 x2  ...  a1n xn  b1 , a x  a x  ...  a x  b ,  21 1 22 2 2n n 2  ... a n1 x1  a n 2 x2  ...  a nn xn  bn приводят последовательным исключением неизвестных к эквивалентной системе с треугольной матрицей  x1  c12 x 2  ...  c1n x n  d1 ,  x 2  ...  c 2 n x n  d 2 ,   ...  xn  d n , решение которой находят по формулам xn  d n , xi  d i  n c k  i 1 ik xk , i  n  1, n  2, ..., 1. Общее число арифметических операций прямого хода метода Гаусса  2 3 n , считая размер3 ность системы n. 0.780 x1  0.563 x 2  0.217 , 0.913 x1  0.659 x 2  0.254 . Пример.  Легко убедиться, что точным решением системы является вектор x  1.000,  1.000 . Если решить эту задачу методом Гаусса, выполняя вычисления с тремя значащими цифрами, то в результате округлений получится совершенно непригодное решение x   0.443, 1.000. Решение системы линейных алгебраических уравнений методом простых итераций Точные методы решения линейных систем применяют для решения линейных систем относительно небольшой размерности (до 10 3 ). Для решения систем большей размерности 3 6 ( 10  10 ) используют итерационные методы. Для решения систем еще больших размерностей(> 106 ) применяют специальные методы, которые здесь не рассматриваются Метод прогонки – алгоритм решения систем линейных алгебраических уравнений с трехдиагональными матрицами. Общее число арифметических операций прямого хода 8n. Итерационные методы хороши для систем с разреженными матрицами. Разреженными называют матрицы, большинство элементов которых — нули. Рассмотрим простейший итерационный метод решения линейной системы — метод простых итераций. Метод состоит в следующем. Сначала система уравнений Cx  d преобразуется к виду x  b  Ax и ее решение вычисляется как предел последовательности x ( k )  b  Ax ( k 1) , k  1, 2, ... Преобразовать систему Cx  d к виду x  b  Ax можно, например, выделив диагональ 1   d i   c ij x j , i  1, 2, ..., n . ные элементы: x i    c ii  i j  Лекция 3. Стр. 6 Теорема. Пусть выполнено условие A  1 . Практически это означает, что итерационные последовательности хорошо сходятся для тех матриц, у которых диагональные элементы велики по сравнению с внедиагональными. Тогда 1) решение x системы x  b  Ax существует и единственно; 2) при произвольном начальном приближении x ( 0) метод простой итерации сходится и справедлива оценка погрешности x (k)  x  A k x ( 0)  x .(*) 1) Из курса линейной алгебры известно, что система линейных алгебраических уравнений имеет единственное решение при любой правой части тогда и только тогда, когда соответствующая однородная система имеет только нулевое решение. Пусть x - решение однородной системы x  Ax . x  Ax  A  x . Т.к. по условию A  1 , это неравенство возможно только при x  0 . Следовательно, x =0 и тем самым первое утверждение теоремы доказано. 2) Вычитая из равенства x ( k )  b  Ax ( k 1) равенство x  b  Ax , получим x ( k )  x  A( x ( k )  x ) . Вычисляя норму левой и правой частей этого равенства и используя неравенство Ax  A  x получим : x (k)  x  A x ( k 1)  x  A 2 x ( k 2)  x  ...  A k 1 x (1)  x  A k x ( 0)  x . Справедливость неравенства (*) установлена. Учитывая, что x ( 0)  x не зависит от k, а A k  0 при k   получаем, что x (k)  x  0 при k   . В качестве условия окончания итерационного процесса можно взять условие x (k )  x ( k 1 ) x e (k )   , где  — заданная погрешность. Тогда полагаем x  x*  x (k) . e Решение системы Cx  d методом Зейделя. В итерационном методе Зейделя (метод Зейделя называют методом Гаусса-Зейделя, процессом Либмана, методом последовательных замещений) xi ( k 1)  1 cii i 1   d i   cij x j ( k 1)   j 1  n c j  i 1 ij xj (k )   . При вычислении каждой последую  щей компоненты решения используются уточненные значения предыдущих компонент. Метод Зейделя легко записать в компактной матричной форме: (k 1) (k 1) (k) x  b  AНx  AВx ,  0   a 21 A   a31   ... a  n1 a12 a13 a32 a 23 ... an2 ... an3 ... a1n   0   ... a 2 n   a 21 0 ... a3n  A   a  н  31 a 32 0 ... ...  ...  ... ...   ... 0   a n1 a n 2 a n 3 Если A Н  A В  1 , то x (k)  x  q k x ( 0 )  x , q  ... 0   0 a12   ... 0  0 0 ,  ... 0 A В   0 0   ... ...  ... ...  0 0 ... 0   AВ 1  A  a13 a 23 ... ... a1n   ... a 2 n  ... a 3n  .  ... ...  ... 0   1 . В практических расчетах Н используют то же условие окончания итерационного процесса, что и в методе простых итераций.