Вычислительная математика лекции. Конспект лекций. Санкт петербург 2011 2 Оглавление
 Скачать 3.86 Mb.
|
|
пограничным слоем. Решение здесь обладает большими производными и изменяется очень быстро. Второй участок продолжительностью много больше первого ( ) nc T обладает сравнительно малыми производными, и решение носит спокойный характер. Для линейных систем 0 0 ( ) ( ), ( ) , x t Ax t x t x с постоянной матрицей свойство "жесткости" предопределяется плохой обусловленностью матрицы 1 ( ) , ( ) 1 cond A A A cond A или max ( ) ( ) , ( ) 1. min ( ) i i i i A k A k A A Для нелинейных систем общего вида 0 0 ( ) ( , ), ( ) x t f t x x t x , жесткость системы часто предопределяется плохой обусловленностью матрицы Якоби f x в процессе решения. На примере решения линейных систем убеждаемся и том, что время наблюдения решения Т определяется минимальным по модулю собственным значением матрицы А, а шаг интегрирования – максимальным. Тогда число шагов N прямо пропорционально числу обусловленности, что приводит к недопустимым затратам max min max min 1 1 , , 1. i i i i T T h N h 123 В связи с указанными причинами следует избегать использования явных методов численного интегрирования. Для интегрирования жестких систем рекомендуется использовать неявный метод ломанных Эйлера или метод Адамса-Мултона второго порядка точности (неявный метод трапеций), устойчивость которых обеспечивается при выполнении условия Re( ( )) 0, i A i , вне зависимости от выбора величины шага интегрирования. Существенный недостаток этих методов малая точность. Среди неявных линейных многошаговых методов присутствуют A(α) устойчивые методы высоких порядков точности, пригодные для решения жестких систем. В качестве примера можно назвать семейство алгоритмов Гира. A(α) устойчивость означает, что область устойчивости метода распространяется на большую часть левой полуплоскости. 11. Решение систем линейных алгебраических уравнений. 11.1. Постановка задачи. Система линейных алгебраических уравнений (СЛАУ), записанная в матричной форме, имеет вид Ax = b, где A R n×n ; x,b R n . Требуется найти её решение с заданной точностью. Прежде чем приступать к решению задачи, необходимо выяснить его принципиальную возможность, т. е. установить существование решения, его единственность и непрерывность зависимости погрешности решения от погрешностей исходных данных (установить корректность задачи ). Затем, рассчитав коэффициенты обусловленности, выяснить возможность 124 получения требуемой точности решения при заданных погрешностях исходных данных (установить практическую возможность решения ). Корректность задачи. Проанализируем условия существования и единственности решения. 1 1 1 n n n x b Ax a a x b Из этой формулы следует, что вектор b есть линейная комбинация столбцов матрицы a i c весами x i (i= 1,2,…,n). Если ранг r(A) =n, т. е. вектор-столбцы линейно не зависимы и, следовательно, матрица не вырождена, то любой конкретный вектор b единственным образом определяется линейной комбинацией вектор столбцов. Иными словами решение системы существует и единственно. В дальнейшем предполагается, что матрица системы не вырождена. В случае вырожденной матрицы r < n, вектор столбцы принадлежат подпространству полного пространства. Это подпространство имеет размерность r. Если вектор b не принадлежит этому подпространству, система несовместна и решение отсутствует. Если вектор b принадлежит подпространству столбцов, то существует бесконечное множество решений. Действительно, общее решение есть сумма решений однородного уравнения и частного решения. Известно, что однородная система с вырожденной матрицей имеет бесконечное число не нулевых решений. Вектор b в этом случае должен иметь n – r элементов, являющихся заданной линейной комбинацией остальных r элементов. Продемонстрируем ситуацию следующим примером. Используя метод исключения Гаусса, получаем 125 1 1 2 2 3 3 1 1 2 3 1 2 2 2 1 2 2 1 0 6 5 0 6 5 1 2 2 0 6 5 2 2 2 1 0 6 5 0 0 0 2 b b b b b b b b b b b b В данном случае r(A)=2. Для существования частного решения необходимо выполнение условия 2b 3 -b 1 -b 2 =0. Неизвестную x 3 полагают равной произвольной константе, x 2 и x 1 определяются в ходе обратной подстановки. Обусловленность задачи. Введем понятия погрешностей вектора и матрицы. Пусть компоненты вектора и матрицы А получают малые приращения. Возмущенная система имеет матрицу А*=А+ΔА и правую часть b*=b+Δb. Ее решение x*=x+Δx. Назовем абсолютной погрешностью вектора и матрицы величины Δ b=|| b*-b ||=|| Δb ||, ΔA=||A* - A||=|| ΔА || соответственно, а величины δb =|| b*-b ||/ || b ||, δA = ||A* - A||/ ||A|| относительными погрешностями. Будет показано, что погрешности решения связаны с погрешностями исходных данных соотношениями || Δx|| = ν Δ (|| Δb || +|| ΔA ||) и || x*-x ||/ || x || = ν δ (|| b*-b ||/ || b || + || A*-A ||/ || A ||). Здесь ν Δ и ν δ коэффициенты абсолютной и относительной обусловленности (чувствительности ) решения задачи. Их величина зависит от свойств матрицы системы и будет вычислена позднее. На заключительном этапе подбирают вычислительный метод по критерию вычислительной устойчивости и эффективности. Процесс решения часто содержит миллионы арифметических операций. Вычислительная неустойчивость приводит к накоплению недопустимо больших погрешностей, а иногда даже к аварийному останову или 126 зацикливанию вычислений. Под эффективностью понимают экономное расходование вычислительных ресурсов. Вычислительные алгоритмы сравнивают по количеству требуемых элементарных операций, по возможности экономного расходования оперативной памяти. Любая хорошо обусловленная задача может быть испорчена вычислительно неустойчивым алгоритмом. С другой стороны самый совершенный алгоритм не справится с плохо обусловленной задачей. При решении СЛАУ используют численные методы двух разновидностей: прямые и итерационные. Характерным признаком прямых методов является возможность получить точное решение (в предположении отсутствия ошибок округления) за конечное число операций. При использовании итерационных методов задают закон, формирующий последовательность векторов x 0 , x 1 , …,x k ,… таким образом, что x k → x при k → ∞, где x точное решение системы. Итерационные методы часто используют при решении больших систем с разреженными матрицами, большинство элементов которых имеют нулевые значения. Такие системы возникают, например, при решении краевых задач в частных производных. Прямые методы, как правило, используют при решении задач среднего размера ( n ≈ 100) с заполненной матрицей. 11.3. Обусловленность задачи решения системы линейных алгебраических уравнений. Предположим сначала, что x решение системы Ax = b и x+Δx решение системы A(x+Δx ) = b+Δb с измененной правой частью. Следовательно, AΔx = Δb, Δx = A -1 Δb, ||Δx|| ||A -1 || ||Δb|| (*). Таким образом, в данных условиях абсолютное число обусловленности ν Δ = ||A -1 ||. С практической точки зрения полезнее использовать относительные погрешности. Из Ax = b следует ||b|| ||A|| ||x||, что вместе с (*) дает 127 |||Δx|| ||b|| A|| ||A -1 || ||Δb|| |||x|| или (b ≠0) (|||Δx|| /|||x||) A|| ||A -1 || (||Δb|| /||b||). Относительное число обусловленности ν δ = ||A|| ||A -1 || . Произведение ||A|| ||A -1 || играет очень важную роль и называется числом обусловленности матрицы A. Стандартное обозначение этого числа cond(A). Рассмотрим теперь случай, когда изменяются элементы матрицы, так что возмущенное уравнение имеет вид (A + ΔA)(x + Δx) = b или с учетом Ax=b AΔx = b - (A + ΔA)x – ΔAΔx =- ΔA(x+Δx ) , -Δx A -1 ΔA(x+Δx ) . Следовательно, || Δx A -1 ΔA (x+Δx )|| = cond(A)( ΔA / A ) ||(x+Δx )||, так что Δx / (x+Δx )|| cond(A)( ΔA / A ) . Обычно ||x|| ≪ Δx||. Поэтому ||(x+Δx )|| ≈ ||(x|| . Если с погрешностью заданы как правая часть, так и матрица, причем cond(A)δ(A) ≪1, можно доказать справедливость неравенства δx ≲cond(A) (δ(b) + δ(A)). Во всех рассмотренных случаях относительное число обусловленности решения совпадает с числом обусловленности матрицы . Полученные результаты нужно правильно интерпретировать. Если cond(A) мало, например, близко к 1, то малые относительные изменения данных обязательно приводят лишь к малому изменению решения. С другой стороны, если число обусловленности велико, то малые изменения в данных могут привести к большим изменениям решения, но это происходит необязательно, а в зависимости от конкретного возмущения. Отметим следующие свойства числа обусловленности. 1. Для единичной матрицы cond(E) = 1. Это следует из того, что E -1 =E и ||E||=1. 128 2. Из условия ||E||=1 и из равенства E=AA -1 следует 1 = ||E|| ||A -1 || ||A|| = cond(A). 3. Число обусловленности матрицы не меняется при умножении матрицы на произвольное число α ≠0, Это следует из соотношений (αA) -1 =α -1 A -1 и ||αA|| = |α| ||A||. На практике влияние большого числа обусловленности зависит от точности данных и длины слова используемой ЭВМ. Проиллюстрируем наши результаты на примере. Рассмотрим систему уравнений с исходными данными 1 1 1.03 0.991 2.51 87.4 91.8 , b= , A 0.991 0.943 2.41 91.8 95.4 ( ) || || || || 2.021*187.2 378 A cond A A A Решением системы является x 1 ≈ 1.981, x 2 ≈ 0.4735. Правая часть системы в лучшем случае известна с точностью 0.005, если считать, что ее элементы получены, например, округлением "истинных" значений при вводе в память трехзначной десятичной ЭВМ. Возмутим каждую из компонент вектора b на величину 0.005, положив b*=(2.505, 2.415) T Решением системы теперь является x*=(2.877, -0.4629) T . Решение оказалось полностью искаженным. Относительная погрешность задания правой части ||b-b*|| ∞ /||b|| ∞ = =0.005/2.51 ≈0.2% привела к относительной погрешности решения ||x –x*|| ∞ / ||x|| ∞ =0.9364/1.981 ≈ 47.3%. Погрешность увеличилась в 237 раз. Это несколько меньше, чем cond(A) = 378. Однако, видимо существует такой вектор b, при котором коэффициент возрастания достигает значения cond(A) . 11.4. Прямые методы решения систем линейных алгебраических уравнений. 11.4.1. Метод Гаусса и LU разложение. 129 1.Схема единственного деления. Из числа прямых наиболее широкое применение находит метод исключения Гаусса. Суть этого метода состоит в эквивалентном преобразовании к системе с верхней треугольной матрицей. При этом формируется последовательность (A | b) →(A 1 | b 1 ) →…→(A k | b k ) →…→ →(A n-1 | b n-1 )=(U | b n-1 ), где U верхняя треугольная матрица. Матрица A k должна содержать нули ниже главной диагонали в первых k столбцах. Пусть уже осуществлено обнуление элементов ниже главной диагонали в предыдущих k-1 столбцах, т. е. сформирована матрица A k-1 . Опишем аналогичную процедуру обработки k – ого столбца, переход (A k-1 | b k-1 ) →(A k | b k ). В простейшем варианте для обнуления элементов k- ого столбца матрицы A k-1 ниже главной диагонали осуществляют следующие операции: рассчитывают коэффициенты l k+i,k = a k+i,k /a k,k (i=1,2,…,n). Здесь a i,j элементы матрицы A k . Строку с k-ым номером матрицы (A k-1 | b k-1 ) умножают на l k+i,k и результат умножения вычитают из строки с i-ым номером этой матрицы (i=1,2,…,n). Указанные операции должны быть последовательно произведены для всех столбцов с номерами k=1,2,…,n-1. В результате (n-1) таких операций получаем систему с верхней треугольной матрицей, последнее уравнение которой содержит лишь одну неизвестную. Последовательной подстановкой снизу вверх (обратный ход) находим остальные неизвестные. Подсчитаем количество необходимых операций. Предположим, что мы договорились назвать каждое деление и каждое умножение -вычитание простой операцией. В самом начале, когда все уравнения имеют n неизвестных, чтобы обнулить n-1 элементов первого столбца ниже главной диагонали, нужно вычислить n-1 штук элементов l и для каждого из них провести (n-1) операций умножения и столько же вычитания, всего Q 1 =2 n(n-1)+n-1=2(n-1) 2 +3(n-1) операций. На каждом следующем шаге 130 число обрабатываемых уравнений уменьшается на единицу. Поэтому Q 2 =2 n(n-2)+n-2=2(n-2) 2 +3(n-2), Q k =2(n-k) 2 +3(n-k) .Общее число операций 1 1 1 2 1 1 1 1 1 2 3 1 1 2 ( ) 3 ( ) 2( 1) (2 1) ( 1) 2 2 3 3 6 2 3 n n n k k k k n n i i Q Q n k n k n n n n n i i n Обратная подстановка идет значительно быстрее. Последнее неизвестное находится с помощью одной операции ( деление на последний ведущий элемент), предпоследнее уравнение требует двух операций (умножение-вычитание, а потом деление). На k–ом шаге требуется k операций, k=1, …, n. Общее число операций 2 1 1 ( 1) 2 2 n k n Q k n n Метод исключения Гаусса эквивалентен разложению матрицы на произведение треугольных матриц A=LU, где U верхняя треугольная матрица, а L нижняя треугольная с единицами на главной диагонали. На первом шаге метода рассчитывают набор коэффициентов l i1 , i=2, …, n. Сформируем матрицу 1 1 21 31 1 0 0 1 0 0 1 L l l . Легко убедится, что 11 12 13 11 12 13 1 1 21 21 22 23 22 21 12 23 21 13 1 31 31 32 33 32 31 12 33 31 13 1 0 0 1 0 0 0 1 0 a a a a a a L A l a a a a l a a l a A l a a a a l a a l a Для обнуления элемента ниже главной диагонали во втором столбце формируем матрицу 1 32 31 12 2 32 22 21 12 32 1 0 0 0 1 0 , 0 1 a l a L l a l a l Тогда 131 11 12 13 1 1 2 1 22 21 12 23 21 13 32 32 31 12 33 31 13 11 12 13 22 21 12 23 21 13 33 31 13 23 21 13 23 1 0 0 0 1 0 0 0 1 0 0 0 0 ( ) a a a L L A a l a a l a U l a l a a l a a a a U a l a a l a a l a a l a l Следовательно, A=LU, где L=L 1 L 2 . Можно показать, что 21 31 32 1 0 0 1 0 1 L l l l Существенным условием возможности такой записи L является нижняя треугольная форма с единичной главной диагональю для матриц 1 1 1 2 и L L Кроме того каждая из этих матриц должна иметь единственный столбец с ненулевыми элементами ниже главной диагонали. Решение системы Ax=b осуществляют в следующей последовательности: 1. Производят факторизацию A=LU. Эта операция требует приблизительно 2n 3 /3 операций. 2. Решают систему с треугольной матрицей Ly=b. Необходимое количество операций приблизительно n 2 /2. 3.Решают систему с треугольной матрицей Ux=y. Необходимое количество операций приблизительно n 2 /2. Общее число операций как и прежде имеет порядок 2n 3 /3. 2.Метод Гаусса с выбором главного элемента по столбцу (схема частичного выбора). В процессе расчета коэффициентов l ij =a ij /a jj , i=j+1, …, n, может оказаться, что ведущий элемент a jj = 0. В этом случае алгоритм становится некорректным, а при условии a jj ≪ a ij вычислительно неустойчивым. Рассмотрим систему 5 1 2 1 10 1 0 2 1 x x с точным решением x 1 =-0.499995…, x 2 =0.999995…. . Пусть длина слова десятичного 132 компьютeра соответствует четырем знакам. Иными словами числа представлены в форме 0.**** ×10 p . Формально мы можем воспользоваться схемой единственного деления, так как a 11 ≠0. Множитель l 21 = 1 6 4 0.2 *10 0.2 *10 0.1*10 вычисляется точно. Значение нового a 22 =0.1*10 1 -(-0.2*10 6 )(0.1*10 1 )=0.1*10 1 +0.2*10 6 =0.2000*10 6 . На самом деле a 22 =0.200001*10 6 , но машинное слово вмещает только четыре цифры. При вычислениях нового b 2 и компонентов x 1 , x 2 в процедуре обратного хода ошибки округления не возникают: b 2 =-(-0.2*10 6 )(0.1*10 1 )=0.2*10 6 6 1 2 6 1 1 1 4 0.2 *10 0.1*10 0.2 *10 0.1*10 0.1*10 0 0.1*10 x x Ошибка округления в шестом знаке после запятой при расчете а 22 привела к полному искажению результата. Отметим, что cond(A)=2, задача хорошо обусловлена. Для исключения подобных ситуаций осуществляют следующую модификацию алгоритма. Во время обработки очередного j-ого столбца просматриваются все его элементы, находящиеся на главной диагонали и ниже ее. Выбирается элемент a ij из условия | | | | max ij ij j i n a a , после чего уравнения с номерами i и j меняются местами. Дальнейшие вычисления идут по прежнему сценарию. Применим модифицированный алгоритм для решения предыдущей задачи. После перестановки строк получаем 1 5 2 2 1 0 10 1 1 x x Решение дает следующие результаты: 4 5 21 1 0.1*10 0.5*10 0.2 *10 l 133 Новое а 21 =0.1*10 1 -(-0.5*10 -5 )(1)=0.1*10 1 (ошибка округления) Новое b 2 =0.1*10 1 -(-0.5*10 -5 )(0)=0.1*10 1 1 1 2 1 1 1 1 0.1*10 0.1*10 0.1*10 (0.1*10 )(1) 0.5 0.2 *10 x x Перестановка i и j строк эквивалентна умножению преобразуемой матрицы слева на матрицу перестановки P ij . P ij получается из единичной матрицы перестановкой ее i и j-ой строк . Для матриц перестановки справедливо соотношение P ij = P ij -1 , P ij P ij =E. Процесс формирования LU разложения модифицируется следующим образом. При обработке очередного k-ого столбца сначала необходимо осуществить перестановку выбранных строк, т. е. сформировать матрицу P k A k . Для этой матрицы найти L k -1 , а затем A k = L k -1 P k A k-1 . Например, цепочка преобразований при n=4 может выглядеть следующим образом A →A 1 →A 2 →A 3 =U. L 3 -1 P 3 L 2 -1 P 2 L 1 -1 P 1 A=U; A 1 =L 1 -1 P 1 A; A 2 =L 2 -1 P 2 A 1 ; A 3 =L 3 -1 P 3 A 2 Так как матрицы P 3 L 2 -1 и P 2 L 1 -1 не являются нижними треугольными с единичной диагональю, получить разложение A=LU из последнего соотношения не удается. Однако это соотношение можно записать в следующей форме L 3 -1 (P 3 L 2 -1 P 3 ) (P 3 P 2 L 1 -1 P 2 P 3 ) ( P 3 P 2 P 1 A)=U . Обозначим 1 1 1 1 1 2 3 2 1 2 3 3 2 3 (P P L P P ), (P L P ) L L , ( P 3 P 2 P 1 A)=PA. Рассмотрим матрицу 1 1 1 3 2 1 2 3 (P P L P P ) L . Здесь L 1 -1 нижняя треугольная матрица с единичной главной диагональю и с единственным, а именно, первым ненулевым столбцом ниже главной диагонали. 134 21 1 1 31 41 1 0 0 0 1 0 0 0 1 0 0 0 1 l L l l Умножение L 1 -1 слева на P 2 преобразует эту матрицу в недиагональную. Вторая строка меняется местами с некоторой, расположенной ниже. Умножение преобразованной матрицы справа на ту же P 2 переставляет столбцы с теми же номерами, что и строки. Результирующая матрица снова оказывается треугольной с единичной главной диагональю, однако со строками, переставленными соответственно структуре P 2 . Аналогичным образом действует преобразование с помощью матрицы P 3 . Таким образом, 1 1 L имеет требуемую форму и отличается от L 1 -1 только перестановкой элементов в первом столбце. То же самое относится к 1 2 L . Разложение имеет вид PA=LU, где L= 1 L 2 L L 3 . Как и прежде не требуется инвертировать и перемножать матрицы. Элементы из столбцов 1 1 L , 1 2 L и 1 3 L просто переносятся с противоположными знаками в те же столбцы матрицы L. Обобщим сказанное на случай матрицы размером n. На первом этапе, используя метод Гаусса с частичным выбором ведущего элемента, получим L n-1 -1 P n-1 L n-2 -1 …P i+1 L i -1 ….P 3 L 2 -1 P 2 L 1 -1 P 1 A=U. В результате определяется набор и последовательность перестановок. Далее, встраивая в нужные места соответствующие сомножители, приводим это разложение к виду 1 1 n L 1 2 n L … 1 i L … 1 1 L P n-1 …P 1 А=U, где 1 1 n L = 1 1 n L , 1 i L = 1 1 1 1 1 n i i i n P P L P P , (i=1,2,…,n-2). 135 Если при обработке очередного k – ого столбца перестановки не требуются, полагаем P k =E. В процессе разложения требуется определить три матрицы P, L и U. Если эти матрицы определены, решение системы разбивается на следующие этапы: 1. Рассчитывают Pb. 2. Решают систему Ly=Pb. 3.Решают систему Ux=y. |