Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
Упражнения. Разработать программу на языке Matlab, реализующую LU разложение без выбора и с частичным выбором ведущего элемента. Применить программу для решения заданной системы линейных алгебраических уравнений. 11.4.2. Метод Холецкого (метод квадратного корня). Метод можно использовать в случае симметрической положительно определенной матрицы А. симметрическая матрица называется положительно определенной, если для любого вектора x выполнено условие x T Ax>0. Симметрическая матрица имеет вещественные собственные числа λ. Она положительно определена, если все её собственные числа λ>0. Линейные системы c матрицами такого типа часто встречаются в приложениях (задачи оптимизации, решение уравнений матфизики). Разложение Холецкого имеет вид 136 11 21 22 1 2 11 21 1 22 2 0 0 0 ; ; 0 0 0 T n n nn n n T nn l l l A LL L l l l l l l l l L l Если разложение выполнено, решение системы Ax=b осуществляется в два этапа: 1. Ly=b; 2. L T x=y. На каждом из этих этапов решается система с треугольной матрицей. Для нахождения элементов матрицы L элементы произведения LL T приравнивают соответствующим элементам матрицы А. Матрица А – симметрическая, поэтому вычисляют только элементы на главной диагонали и ниже её ( ; , ; 1, ) ij a i j n j n . В результате получают ( 1) 2 n n уравнений для того же количества неизвестных l ij . 2 11 11 1 11 1 2 2 21 22 22 1 21 2 22 2 2 2 2 1 2 1 1 2 2 ; i=2, ; элемeнты 1 -ого столбца ; i=3, ; элемeнты 2 -ого столбца i i i i i k k kk kk i k i k ik l a l l a n l l a l l l l a n l l l a l l l l l 2 2 2 1 2 ; i= 1, элемeнты k -ого столбца ; элемeнты n -ого столбца kk ik n n nn nn l a k n l l l a Эта система решается последовательно: сначала вычисляют l 11 затем l i1 , l 22 и т. д. 137 11 11 1 1 11 2 22 22 21 2 2 1 21 22 2 2 2 1 2 , 1 1 1 2 2 , 1 , 1 / ; i=2, ( ) / ; i=3, ; 1, ( ) / ; i= 1, i i i i i kk kk k k k k ik ik i k i k i k k k kk l a l a l n l a l l a l l l n l a l l l k n l a l l l l l l l k n 2 2 2 1 2 , 1 nn nn n n n n l a l l l Подкоренные выражения положительны, если матрица А положительно определена. Компактная форма записи тех же расчетных формул: 1. Для 1, k n проделать 1 2 1 1 1 1.1 ; 1.2 Для i= 1, проделать k kk kk ki i k ik ij kj j ik kk l a l k n a l l l l 2. Вывод матрицы L. Необходимое число операций Q ≈n 3 /3, что вдвое меньше требуемого в методе Гаусса. Другим положительным качеством метода является высокая вычислительная устойчивость алгоритма. Упражнения. Разработать программу на языке Matlab, реализующую метода Холецкого . Применить программу для решения заданной системы линейных алгебраических уравнений. 138 11.4.3.QR разложение матрицы и решение систем линейных алгебраических уравнений. QR разложение означает факторизацию вида A=QR, где Q – ортогональная матрица, R – верхняя треугольная матрица. Матрица Q является ортогональной, если выполнено хотя бы одно ( а, следовательно, все) из следующих условий: 1. Q -1 = Q T ; 2. Cтолбцы Q образуют ортонормированную систему (т. е. (q i , q j )= 0 и (q i , q i ) = 1); 3. Cтроки Q, рассматриваемые как векторы, образуют ортонормированную систему; 4. ||Qx|| = ||x|| для любого x. Из второго и третьего условия следует, что если Q ортогональная матрица, то и Q T ортогональная. Из условия (4) следует, что ||Q|| p =1, а из (1) и (4) – cond(Q)=1 (cond(A)=||A|| ||A -1 ||). Отметим следующее свойство: произведение ортогональных матриц есть матрица ортогональная. Если вычислено A =QR, система Ax = b решается следующим образом. 1. Решают систему Qy = b, y =Q T b. При этом отпадает необходимость обращения матрицы Q. 2. Решают систему с верхней треугольной матрицей Rx = y. 1. Метод вращения. Матрица вращения T kl получается из единичной матрицы добавлением четырех элементов, расположенных на пересечении k-ой и l-ой строк с k-ым и l-ым столбцами (k 139 1 1 1 kl kl kl kl kl c s T s c Все элементы, кроме отмеченных, равны нулю. Матрица вращения ортогональна, т. е.c kl 2 +s kl 2 =1, что эквивалентно s kl =sin(φ kl ), c kl =cos(φ kl ). Вектор y=T kl x отличается от вектора x элементами в k-ой и l-ой строках y k =c kl x k +s kl x l , y l =c kl x l -s kl x k . Действие T kl на вектор x эквивалентно его повороту на угол φ kl вокруг оси, перпендикулярной гиперплоскости Оx k x l . Для обнуления элемента x l требуется выполнение следующих условий: y l =c kl x l -s kl x k =0, c kl 2 +s kl 2 =1, откуда следует 2 2 2 2 , s k l kl kl k l k l x x c x x x x Преобразование исходной матрицы A к верхней треугольной R осуществляется в результате последовательного обнуления элементов столбцов, расположенных ниже главной диагонали, начиная с первого столбца и заканчивая (n-1)-ым: A →A 1 →…→A k →…→A n-1 =R . Первые k столбцов матрицы A k имеют нулевые компоненты ниже главной диагонали. На первом шаге вычисляют A 1 , обнуляя (n-1) элементов a i1 (i=2,…, n) первого столбца матрицы А. Для исключения элемента a 21 формируют матрицу вращения T 12 , полагая x k =a 11 , x l =a 21 , и вычисляют T 12 A=(T 12 a 1 ,T 12 a 2 ,…T 12 a i ,…,T 12 a n ). Здесь a i i-ый столбец матрицы А. Для исключения очередного элемента a i1 первого столбца необходимо определить матрицу вращения T 1i , положив x k =a 11 , x l =a i1 . В результате получим A 1 =T 1n T 1,n-1 …T 13 T 12 A=T 1 A. Матрица T 1 ортогональна, так как является произведением ортогональных матриц. На k-ом шаге при 140 переходе A k-1 →A k после обнуления элементов k-ого столбца получают A k = T k,n T k,n-1 …T k,k+2 T k,k+1 A k-1 =T k A k-1 . Выполнив (n-1) шагов (k=1,2,…,n-1), получим T n-1 T n-2 …T 1 A=R или A=QR, где Q= T 1 T T 2 T …T n-1 T При формировании соответствующих матриц вращения учитывают номер обнуляемого элемента. Если, например, обнуляется элемент a ij (i>j), то элементы матрицы T ji равны (x k =a jj , x l =a ij ). 2 2 2 2 , s jj ij ji ji jj ij jj ij a a c a a a a Количество требуемых операций Q ≈2n 3 . Метод обладает повышенной вычислительной устойчивостью в сравнении с методом Гаусса. 2.Метод отражения. Матрица отражения (матрица Хаусхольдера) V=E-2ww T , w R n , V R n×n , ||w||= 2 2 2 1 2 T n w w w w w =1. V= 2 1 1 2 1 2 2 1 2 2 2 1 2 1 2 2 2 2 1 2 2 2 2 1 2 n n n n n w w w w w w w w w w w w w w w Матрица отражения симметрическая и ортогональная. Действительно V=V T и V T V=(E-2ww T )(E-2ww T )=E-4ww T +4ww T ww T =E. То есть V T =V -1 В основе использования матрицы отражения для QR разложения лежит следующая теорема Теорема. Пусть заданы любые ненулевые векторы a и b. Существует единственный вектор w единичной длины такой, что определяемая им матрица отражения переводит вектор a в b (Va=b). Действие матрицы V на вектор a есть зеркальное отражение вектора a относительно гиперплоскости, проходящей через начало координат и перпендикулярной вектору w. 141 A 1 =V 1 A=V 1 (a 1 a 1 … a 1 )=(V 1 a 1 V 1 a 2 … V 1 a n ). В соответствии с теоремой всегда можно построить V 1 таким образом, чтобы V 1 a 1 имел нули ниже главной диагонали. Далее подбираем V 2 так, чтобы во втором столбце V 2 A 1 обнулить элементы ниже главной диагонали. В этих условиях будут сохранены нули ниже главной диагонали и в первом столбце. Проделав эти операции (n-1) раз, получим V n-1 V n-2 …V 1 A=R . V n-1 V n-2 …V 1 =Q T ортогональная матрица. Q T A=R или A=QR. Правило выбора вектора w k . (V k =E-2w k w T k ). 1, 2, , 1 2 0,0,...0, , , ,..., ; 1 ; ; k=1, 1. 2 ( ) T k k kk k k k k k n k k нулей n k ik k i k k k kk w a s a a a s a n s s a Знак s k совпадает со знаком –a kk Количество требуемых операций Q ≈ 4/3)n 3 . Метод обладает повышенной вычислительной устойчивостью в сравнении с методом Гаусса. Упражнения. Разработать программу на языке Matlab, реализующую QR разложение с использованием метода отражения. Применить программу для решения заданной системы линейных алгебраических уравнений. 11.5. Сравнительный анализ вычислительной устойчивости прямых методов решения систем линейных алгебраических уравнений. 142 Метод Число операций Точность (f(n)) Доп. память Примечания Метод Гаусса 2/3n 3 αn 2 0 α< ∞(ед. деление) α< 2 n-1 (част. выбор). реально α< ÷ n Метод Холецкого 1/3n 3 1 0 Метод вращения 2n 3 2.9n 0 Метод отражения 4/3n 3 2.9n 2n Вычислительная погрешность оценивается по формуле || * || (2 ) ( ) || || M x x condA f n x 11.6. Итерационные методы решения систем линейных алгебраических уравнений. Итерационный метод задает процедуру построения последовательности векторов x 1 , x 2 , …, x n , … такую, что x n →x* при n→∞, где x* точное решение исходной системы Ax=b. Построение последовательности прерывается при достижении заданной погрешности ε ||x n -x*||<ε. Указанная процедура может быть сформирована следующим образом. Эквивалентным преобразованием исходная система приводится к виду x=Cx+d , (1) после чего заменяется разностной системой 143 x n+1 =Cx n +d (2). Требуют ответа следующие вопросы. 1. Как осуществить переход к системе (1). 2. Сходится ли итерационный процесс, а если сходится, то к какому пределу. 3.Какова скорость сходимости. Начнем с ответа на второй вопрос. Точное решение исходной системы удовлетворяет уравнению (1) x*=Cx*+d. Вычитая последнее уравнение из (2), получаем e n+1 =Ce n (3) e n =x n -x* погрешность на n-ом шаге. Для того, чтобы x n →x* при n →∞, требуется e n →0 при n→∞. Система (3) должна быть асимптотически устойчивой. Необходимое и достаточное условие сходимости итерационного процесса | ( ) | 1. max i i C Из за сравнительной сложности вычисления собственных чисел предпочитают использовать достаточное условие сходимости ||C||<1. Сходимость в этом случае доказывается следующим образом. Из (3) следует e n =C n e 0 , ||e n || C|| n ||e 0 ||. Таким образом, при n →∞ ||e n || →0, если ||C||<1. Для проверки сходимости можно использовать любую норму матрицы. Количество необходимых итераций тем меньше (тем больше скорость сходимости), чем меньше модули собственных чисел матрицы С и чем ближе к точному решению выбрано начальное приближение x 0 Обсудим условия окончания итерационного процесса. Теорема. Если выполнено условие ||C|| <1, то справедлива следующая оценка погрешности 1 || || || *|| || || . 1 || || n n n C x x x x C Доказательство. Имеем x n -x*=C(x n-1 -x*). 144 Иначе x n -x*=C(x n-1 -x*)+C(x n -x*). Тогда || x n -x*||=||C(x n-1 -x*)+C(x n -x*)|| ||C|| ||(x n-1 -x*)||+||C|| ||(x n -x*)|. Таким образом, 1 || || || *|| || || . 1 || || n n n C x x x x C При переходе от исходной системы к системе (1) требуется сформировать матрицу С с указанными свойствами. Универсальный способ такого перехода отсутствует. В каждом конкретном случае требуется анализировать индивидуальные особенности исходной системы. Рассмотрим некоторые из широко используемых приемов формирования системы (1). Метод Якоби. Пусть диагональные элементы матрицы А отличны от нуля. Тогда исходную систему можно представить в виде ( ) 1 ( ) ( ) ( ) 1 1 , 1, 2,..., . k k m kj kj k j j j j k kk kk kk a a b x x x k m a a a Здесь А R m×m , x (k) – k-ый компонент вектора x. Переходя к разностной системе, получаем ( ) 1 ( ) ( ) ( ) 1 1 1 , 1, 2,..., . k k m kj kj k j j n n n j j k kk kk kk a a b x x x k m a a a (4) Это метод Якоби. В матричной форме получаем x n+1 =-D -1 (A-D)x n +D -1 b или D(x n+1 -x n )+Ax n =b. Здесь D диагональная матрица с диагональными элементами матрицы А, C=-D -1 (A-D). Чем больше абсолютные величины диагональных элементов матрицы А по сравнению с модулями остальных ее элементов, тем меньше модули элементов матрицы С и ее норма и тем выше скорость сходимости итерационного процесса. Метод Зейделя. К моменту вычисления x n+1 (k) значения x n+1 (j) для j 145 ( ) 1 ( ) ( ) ( ) 1 1 1 1 , 1, 2,..., . k k m kj kj k j j n n n j j k kk kk kk a a b x x x k m a a a Это метод Зейделя. Введем диагональную матрицу D, а также нижнюю и верхнюю треугольные матрицы А 1 и А 2 с нулевыми главными диагоналями. Элементы этих матриц совпадают с соответствующими элементами матрицы А, т. е. A=D+A 1 +A 2 . матричная форма метода Зейделя имеет вид x n+1 =-D -1 A 1 x n+1 - D -1 A 2 x n +D -1 b или (D+A 1 )(x n+1 -x n )+Ax n =b. Для метода Зейделя C=(D+A 1 О сходимости методов Якоби и Гаусса – Зейделя. |