Методы решения систем линейных уравнений (PDF)

3 Оглавление ВВЕДЕНИЕ.................................................................................................. 5 1. СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ И МЕТОДЫ ИХ РЕШЕНИЯ......................................................................................

3 Оглавление ВВЕДЕНИЕ.................................................................................................. 5 1. СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ И МЕТОДЫ ИХ РЕШЕНИЯ.................................................................................................... 6 1.1. Метод Гаусса........................................................................... 11 1.2. Метод простых итераций....................................................... 15 1.3. Сравнительная оценка прямых и итерационных методов.. 20 2. НАХОЖДЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИЙ И СОБСТВЕННЫХ ВЕКТОРОВ МАТРИЦ................................................ 22 2.1. Определение наибольшего и наименьшего собственных значений итерационным методом....................................................... 25 3. АППРОКСИМАЦИЯ И ИНТЕРПОЛЯЦИЯ ФУНКЦИЙ.................. 30 3.1. Постановка задачи и основные определения............................... 30 3.2. Интерполяция с помощью многочленов...................................... 34 3.3. Интерполяционный многочлен Лагранжа................................... 34 3.4. Интерполяционный многочлен Ньютона.................................... 37 3.5. Точность и сходимость многочленной интерполяции............... 39 3.6. Использование локальных интерполяций................................... 41 4. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ.......................................... 43 4.1. Постановка задачи......................................................................... 43 4.2. Использование ряда Тейлора........................................................ 45 5. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ................................................... 48 5.1. Постановка задачи......................................................................... 48 5.2. Метод прямоугольников............................................................... 51 5.3. Метод трапеций.............................................................................. 53 5.4. Метод Симпсона............................................................................ 56 6.ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ И СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ................. 59 6.1. Нелинейные уравнения................................................................. 59 6.1.1. Постановка задачи.................................................................. 59 6.1.2. Отделение корней нелинейного уравнения......................... 61 6.1.3. Метод половинного деления и метод хорд.......................... 62 6.1.4. Метод простой итерации....................................................... 66 6.1.5. Метод Ньютона...................................................................... 69 6.1.6. Метод секущих....................................................................... 71 6.1.7. Сравнительная оценка методов............................................. 72 6.2. СИСТЕМЫ НЕЛИНЕЙНЫХ УРАВНЕНИЙ............................... 74 6.2.1. Постановка задачи и ее особенности.................................... 74 6.2.2. Метод Ньютона и его модификации.................................... 77  4 7. ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ......... 83 7.1. Общие понятия............................................................................... 83 7.2. Метод конечных разностей........................................................... 86 7.3. Решение задачи Коши.................................................................... 89 7.4. Метод Эйлера................................................................................. 91 7.5. Метод Эйлера с пересчетом.......................................................... 92 7.6. Метод Рунге–Кутта........................................................................ 94 7.7. Многошаговые методы.................................................................. 97 Лабораторная работа № 1 Решение систем линейных алгебраических уравнений точными методами.................................... 100 Лабораторная работа № 2 Решение систем линейных алгебраических уравнений численными методами.......................................................... 103 Лабораторная работа № 3 Вычисление собственных значений и векторов матриц....................................................................................... 105 Лабораторная работа № 4 Интерполирование функций...................... 108 Лабораторная работа № 5 Дифференцирование функций................... 110 Лабораторная работа № 6 Приближенные методы вычисления собственных интегралов......................................................................... 112 Лабораторная работа № 7 Аппроксимация функций методом наименьших квадратов............................................................................ 118 Лабораторная работа № 8 Приближенные методы решения нелинейных уравнений и систем нелинейных уравнений......................................... 121 Лабораторная работа № 9 Численные методы решения обыкновенных дифференциальных уравнений............................................................... 125 Заключение............................................................................................... 128 Библиографический список.................................................................... 128  5 ВВЕДЕНИЕ В практической деятельности часто приходится иметь дело с ма- тематическими задачами, точное решение которых является довольно сложным или неизвестно. В этих случаях выходом из создавшейся си- туации являются приближенные вычисления. По этой причине при- ближенные методы решений имеют большое значение и широкое при- менение. В настоящее время эффективное использование ЭВМ для решения научных и инженерных задач практически невозможно без достаточ- ных знаний вычислительной математики. Кроме того, следует отме- тить, что при приближенном решении задач возникает вопрос необхо- димости оценки погрешности полученного результата. Для практического усвоения рассматриваемого в курсе лекций тео- ретического материала авторами настоящего издания составлены соот- ветствующие лабораторные работы. К каждой лабораторной работе приводится необходимый теоретический материал, обращение к кото- рому позволит успешно выполнить индивидуальное задание.  6 1. СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ И МЕТОДЫ ИХ РЕШЕНИЯ Решение систем линейных алгебраических уравнений (СЛАУ) и нахождение собственных значений матрицы относят к основным зада- чам линейной алгебры. Численные методы, применяемые при реше- нии вышеперечисленных задач, являются основой вычислительной математики в силу того, что: – эти задачи достаточно часто фигурируют в научных, инженерных и экономических расчетах; – к ним сводится решение более сложных математических задач; – при построении численных методов решения нелинейных задач кон- струируются алгоритмы, базирующиеся на решении линейных задач. Система n линейных алгебраических уравнений с n неизвестными в общем виде может быть записана следующим образом a11x1  a12 x2 ...  a1n xn  b1 ,  a21x1  a22 x2 ...  a2 n xn  b2 , (1.1)  ................................................. an1 x1  an 2 x2 ...  ann xn  bn. Систему (1.1) часто записывают в компактном матричном виде   Ax  b , (1.2) где использованы следующие обозначения:  a11a12...a1n   x1   b1         a a...a  ;   x  ;   b2 . A   21 22 2 n  x   2  b   .....................        a a...a  x  b   n1 n 2 nn   n  n  Здесь A – матрица коэффициентов при неизвестных системы; b – вектор-столбец заданных правых частей; x – вектор-столбец искомого решения. Система (1.1) называется однородной, если все числа, стоящие в правой части системы, равны нулю (bi = 0, i = 1, 2, … , n). В противном случае, если существует хотя бы одно bi отличное от нуля, система называется неоднородной. Геометрический смысл произведения мат-  рицы на вектор A заключается в том, что, подействовав на вектор x   матрицей A , мы получаем другой вектор y  Ax , который будет по-  вернут на некоторый угол относительно x и длина которого изменится  по сравнению с длиной вектора x. Таким образом, геометрически ре-  7 шение системы (1.2) заключается в том, чтобы найти такой вектор x ,  который матрица A переводит в вектор b. Различные методы решения системы (1.1) основываются на пре- образовании ее в какой-либо другой вид, откуда решение может быть легко получено. Конечно, преобразования при этом должны быть та- кими, чтобы они не изменяли решения системы. Преобразования ли- нейной системы уравнений, не изменяющие решение системы, назы- ваются эквивалентными преобразованиями. Общий вид таких преобра- зований заключается в том, что обе части уравнения системы в мат- ричном виде (1.2) умножаются на одну и ту же невырожденную мат- рицу (определитель которой не равен нулю)   BAx  Bb.   Вводя матрицу C = BA и новый вектор q  Bb , получим исходную систему (1.2) в преобразованном эквивалентном виде:   Cx  q. (1.3) Преобразования имеют смысл, если полученную систему решить гораздо легче, чем исходную систему (1.2). Во многих случаях наряду с матрицей A системы удобно рассматривать расширенную матрицу, которая получается добавлением к основной матрице вектора правых частей системы  a11a12...a1nb1     a21a22...a2 nb2  (1.4) A . .....................   a a...a b   n1 n 2 nn n  В теории линейной алгебры формулируются условия существова- ния решения для линейной системы (1.1). Система (1.1) имеет решение тогда и только тогда, когда ранг матрицы A и ранг расширенной мат- рицы (1.4) совпадают. Если при этом ранг матрицы A равен n (это означает, что определитель матрицы A не равен нулю и система явля- ется неоднородной), то система (1.1) имеет единственное решение, определяемое по правилу Крамера: xi  Di / D , (1.5) где D – определитель матрицы A, а Di – определитель матрицы, полу- ченной из A заменой i-го столбца на столбец правых частей b. Однако простота формулировки способа нахождения решения является обман- чивой. Все существенно осложняется при переходе к численному ре- шению линейных систем на ЭВМ. Применение правила Крамера на ЭВМ при больших n чрезвычайно затруднено в связи со следующими причинами:  8 – наблюдается большой рост вычислительной сложности алгоритма, например, для вычисления одного определителя путем вычисления всех входящих в сумму членов требуется около n  n! арифметических операций; уже при n >10 время его работы становится практически нереализуемым; – увеличение числа операций приводит к быстрому росту погрешности округления и может привести к переполнению разрядной сетки ком- пьютера. Кроме того, на практике приходится решать системы достаточно высокого порядка (от нескольких сотен до сотен тысяч уравнений), что приводит к недопустимому росту погрешности, если не использовать специальных приемов в процессе решения. Если ввести понятие обратной матрицы на основании соотноше- ния 1 A A E, (1.6) где A – исходная матрица; A–1 – обратная к исходной; E – единичная матрица, можно легко получить решение линейной системы (1.2). Действительно, умножив обе части (1.2) слева на A–1, получим 1  1   1   1  A Ax  A b  Ex  A b  x  A b. Следовательно, умножив матрицу A–1 на вектор правых частей b, мы получаем решение системы (1.2). Вообще говоря, нахождение об- ратной матрицы – задача более сложная, чем решение системы (1.2). Находить ее имеет смысл в случае, когда она представляет самостоя- тельный интерес для решаемой задачи, либо при необходимости нахо- дить ряд решений системы (1.2) с различными правыми частями. Методы решения систем линейных уравнений можно разделить на две группы – прямые (или точные) и итерационные. Прямые методы – это методы, которые при отсутствии погрешностей округле- ния за конечное число арифметических и логических операций дают точный результат. Эти методы сравнительно просты. Их иногда назы- вают точными, поскольку решение выражается в виде точных формул через коэффициенты системы и правые части. Однако точным решение может быть лишь теоретически, при выполнении вычислений с беско- нечным числом разрядов. На практике при использовании ЭВМ вы- числения проводятся с ограниченным числом разрядов и поэтому неизбежны погрешности в окончательных результатах решения задачи. Итерационные методы заключаются в построении последова- тельности приближенных решений, предел которой (если он существу- ет) является решением системы.  9 Выбор конкретного численного метода определяется специфиче- скими особенностями решаемой задачи и располагаемой вычислитель- ной техникой. В первую очередь здесь необходимо учитывать вид мат- рицы решаемой линейной системы. Приведем наиболее типичные ви- ды матриц, которые встречаются в различных практических задачах и являются основой многих методов. Диагональная матрица – квадратная матрица, у которой aij = 0 при i  j (i = 1, 2, … n; j = 1, 2, … n). Единичная матрица – диагональная матрица, у которой на глав- ной диагонали стоят единицы, т.е. aii = 1 (i = =1, 2, … n), aij = 0 при i  j. Верхняя треугольная матрица – матрица, у которой ниже глав- ной диагонали все элементы равны нулю, т.е. aij = 0 при i > j. Нижняя треугольная матрица – матрица, у которой выше глав- ной диагонали все элементы равны нулю, т.е. aij = 0 при i < j. Симметричная матрица – квадратная матрица, у которой все элементы, симметричные относительно главной диагонали, равны, т.е. aij = aji (i = 1, 2, …n; j = 1, 2, … n). Ленточная матрица – квадратная матрица, у которой все ненуле- вые элементы образуют “ленту”, параллельную главной диагонали, а все остальные элементы матрицы равны нулю, т.е. aij =0 при i  j  l , где l – ширина ленты. Трехдиагональная матрица – квадратная матрица, у которой не равны нулю только те элементы, которые расположены на главной диагонали и двух соседних с ней диагоналях, т.е. ai j 0 при i  j  1. Трехдиагональная матрица – частный случай ленточной матрицы. Разреженная матрица – квадратная матрица, у которых количе- ство ненулевых элементов гораздо меньше общего количества элемен- тов. В табл. 1.1 приведены примеры описанных матриц. Если рассмотреть матрицу как упорядоченный набор чисел, кото- рый нужно хранить в памяти компьютера, то их можно разделить на два класса – хранимые матрицы и воспроизводимые матрицы. Матрица называется хранимой, если при проведении расчетов ее элементы необходимо хранить в памяти компьютера. Матрица называ- ется воспроизводимой (вычисляемой), если ее элементы не требуется хранить в памяти, так как в процессе вычислений по мере необходимо- сти элементы матрицы могут быть вычислены по простым формулам (воспроизведены). Примером воспроизводимой матрицы является сле- дующая трехдиагональная матрица:  10 5 1 0 0 0 –1 5 1 0 0 А= 0 –1 5 1 0 0 0 –1 5 1 0 0 0 –1 5 Ее элементы легко вычисляются: aii = 4 на главной диагонали, т.е. при i = j; aij = 1 при j = i + 1 (выше главной диагонали); aij = −1 при j = i – 1 (ниже главной диагонали); aij = 0 в остальных случаях. Таблица 1.1 Примеры различных матриц Тип матрицы Пример Диагональная 5 0 0 0 0 −7 0 0 0 0 8 0 0 0 0 4 Единичная 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 Верхняя треугольная матрица 5 1 0.8 5 0 −7 1 2 0 0 8 0 0 0 4 Нижняя треугольная матрица 5 0 0 0 −7 0 0 3 4 8 0 1 9 −4 4 Симметричная матрица 5 10 5 7 10 - 0 9 5 0 8 0 7 9 0 4 Ленточная матрица 5 1 2 0 0 0 9 −7 5 1 0 6 7 8 1 2 0 0 4 5 4 3 0 0 0 4 5 6 1 0 0 0 5 8 4  11 Окончание табл. 1.1 Тип матрицы Пример Трехдиагональная матрица 5 1 0 0 0 9 −7 5 0 0 0 7 8 1 0 0 0 5 4 3 0 0 0 1 6 Разреженная матрица 5 0 0 0 0 0 0 −7 0 1 0 0 0 0 0 0 2 0 0 0 5 0 3 0 0 0 0 0 0 1 0 0 0 0 8 4 1.1. Метод Гаусса Метод Гаусса – это один из наиболее эффективных прямых мето- дов решения систем линейных уравнений. Согласно этому методу ис- ходную систему (1.2) с помощью эквивалентных преобразований необходимо привести к треугольному виду  * Ux  b , (1.7) * где U – верхняя треугольная матрица, а b – вектор правых частей с учетом изменений в процессе эквивалентных преобразований. По- дробнее (1.7) выглядит следующим образом:  u11x1  u12 x2 ...  u1,n1 xn1  u1n xn  b1 , *   u22 x2 ...  u2,n1 xn1  u2 n xn  b2* , (1.8)  ....................................................................   un1,n1 xn1  un1,n xn  bn*1 ,   unn xn  bn*.  Система (1.8) удобна тем, что из нее можно непосредственно опре- делить решение с помощью простого алгоритма. Действительно, из последнего уравнения сразу находим xn. Подставив полученное значе- ние в предпоследнее уравнение, мы получим уравнение с одним неиз- вестным xn-1, откуда сразу определяется его значение. Далее последо- вательно определяются значения искомого решения xn-2, xn-3, …, x2, x1.  12 Последним будет найдено значение неизвестной x1 из первого уравне- ния системы (1.8). Таким образом, метод Гаусса состоит из двух шагов: прямого и об- ратного хода. Прямой ход заключается в сведении системы (1.2) к тре- угольному виду (1.8). Обратный ход заключается в последовательном вычислении искомых переменных из системы вида (1.8). Пример 1.1. Методом Гаусса решить систему уравнений 2 x1  x2  x3  3,  (1.9)  x1  x2  2 x3  1, 3x  2 x  x  8.  1 2 3 Для того чтобы исключить из второго и третьего уравнений пере- менную х1 , из второго уравнения почленно вычтем первое уравнение, умноженное на 1/2, а из третьего уравнения вычтем первое, умножен- ное на 3/2. Такие преобразования приводят к эквивалентной системе, имеющей то же решение, что и исходная система. Подчеркнем, что при умножении уравнения на число производится умножение как правой, так и левой части уравнения. Преобразованная система имеет вид  2 x1  x2  x3  3,  (1.10)  3 5 1 0  x2  x3   ,  2 2 2  1 5 7 0  2 x2  2 x3  2. Для того чтобы привести систему к треугольному виду, нам оста- лось исключить х 2 из третьего уравнения (1.10). Для этого из третьего уравнения вычтем второе уравнение, умноженное на −1/3. В результа- те получим  2 x1  x2  x3  3,  (1.11)  3 5 1 0  x2  x3   ,  2 2 2  10 10 0  0  3 x3  3. На этом заканчивается прямой ход метода Гаусса, так как матрица системы (1.11) является верхней треугольной. С помощью обратного хода находим решение: из третьего уравнения находим х3  1. Подста- вив это значение во второе уравнение, найдем, что х2  2. Подстав-  13 ляя х3  1 и х2  2 в первое уравнение, определим, что х1  1. Итак, решением системы уравнений (1.9) является вектор с координатами х1  1 , х2  2 , х3  1. Рассмотрим теперь прямой ход в общем случае для системы (1.1). Во всех уравнениях, кроме первого, исключается переменная x1. Для этого из второго уравнения вычтем первое, умноженное на а21 а11 , из третьего – первое, умноженное на а31 а11 и так далее, т. е. из i-го уравнения вычтем первое, умноженное на аi1 а11. В результате исход- ная система сведется к виду а11 x1  а12 x2 ...  а1n xn  b1 ,   (1) (1) а22 x2 ...а2 n xn  b2 , (1)  (1.12)  .................................................  (1) (1) (1)  an 2 x2 ...  ann xn  bn. Верхний индекс коэффициентов и правых частей во всех уравне- ниях, начиная со второго, указывает на то, что эти значения коэффици- ентов новые, полученные в результате преобразования системы (1.1). При преобразовании системы мы предполагали, что a11  0. На следу- ющем шаге исключим переменную х2 из всех уравнений, кроме перво- го. Для этого, выбирая в качестве ведущего элемента а 22 (1) , вычтем из третьего уравнения (1.12) второе уравнение, умноженное на а32 (1) а22 (1) и так далее, т. е. из i-го уравнения вычтем второе уравнение, умножен- ное на аi(21) а22 (1) для всех 3  i  n. После таких преобразований систе- ма примет вид а11x1  а12 x2  a13 x3 ...  а1n xn  b1 ,   а22 (1) x2  a23 (1) x3 ...  а2(1n) xn  b2(1) ,  (1.13)   a33( 2) x3 ...  а3( 2n) xn  b3( 2 ) ,  ........................................................   an( 23) x3 ...  ann( 2) xn  bn( 2 ).   14 Описанный процесс продолжается до тех пор, пока исходная си- стема (1.1) не примет треугольный вид (1.8). При описании прямого хода метода Гаусса предполагалось, что на каждом шаге в качестве ведущих были приняты элементы, стоящие на главной диагонали и не равные нулю. Покажем, что это ограничение не существенно. Действительно, если, например, a11 = 0, то мы выберем такое i-е уравнение, в котором коэффициент при x1 не равен нулю ( аi1  0 ). Та- кое уравнение обязательно найдется, так как в противном случае весь первый столбец матрицы системы (1.1) состоит из нулей и, следова- тельно, определитель этой системы равен нулю и матрица A является вырожденной. Поменяем теперь местами первое и i-е уравнения и, приняв за ведущий элемент ai1, проделаем все описанные выше и исключим x1 из всех уравнений, кроме i-го (которое теперь является первым). Аналогично можно поступать и в случае, когда какой-либо ведущий элемент на другом шаге равен нулю. Существуют различные видоизменения метода Гаусса, которые, сохраняя его основную идею, отличаются способом хранения матриц коэффициентов в памяти ЭВМ, способом исключения элементов, ме- тодом ”борьбы” с погрешностями округления. Имеются различные модификации метода для решения задач с матрицами специального вида (симметричными, ленточными и т.д.). Доказано, что метод Гаусса – один из наиболее экономичных пря- мых методов для систем общего вида с точки зрения количества требу- емых операций. В общем случае для решения системы n уравнений требуется приблизительно n 3 / 3 операций. Для понимания возможно- стей метода Гаусса оценим на примере влияние порядка системы n на общее время решения системы и величину необходимой оперативной памяти. Посчитаем, сколько времени понадобится для решения 104 уравнений с плотно заполненной матрицей на компьютере со средним быстродействием 106 операций в секунду. Требуемое количество арифметических операций приблизительно будет равно N  1 3 1012 , следовательно, время решения задачи на компьютере примерно равно t  1 3 1012 / 106 с  1 3 106 с  93 ч  3,8 сут. Для хранения матрицы потребуется примерно 4  n 2 байтов, что в рассматриваемом случае составляет 4108 байт=400 гигабайт (предпо- лагается, что для хранения одного числа используется 4 байта). Для сравнения укажем, что объем памяти винчестера компьютера Pentium2 может варьироваться от 2 до 10 гигабайт, следовательно, решение по- добной задачи на таких компьютерах не представляется реальным. Для  15 случая линейной системы при n = 103 мы получим время счета поряд- 9 ка t  1  106 с  0,1 ч  6 мин и объем памяти 4106 байт = 4 гигабайта. 3 10 Эти оценки показывают, что на современных персональных компьюте- рах методом Гаусса можно решать линейные системы с плотно запол- ненной матрицей порядка n1000. С другой стороны, во многих науч- ных и инженерных задачах приходится решать системы порядка десят- ков и сотен тысяч уравнений. Это удается делать потому, что эти мат- рицы в подавляющем большинстве случаев являются сильно разре- женными и воспроизводимыми, что позволяет для их решения исполь- зовать либо специальные модификации метода Гаусса либо итераци- онные методы. Один из наиболее часто применяемых – метод Гаусса с частичным выбором ведущего элемента, идея которого состоит в том, что в качестве ведущего элемента в каждом рассматриваемом столбце выбирается максимальный по модулю элемент в этом столбце. Опи- санный выше алгоритм прямого хода метода Гаусса фактически реали- зуется следующим образом. На первом шаге в первом столбце отыски- вается максимальный по модулю элемент. Пусть, например, это будет элемент ai1. Поменяем i-е и первое уравнения местами, так что макси- мальный элемент окажется в первой строке и первом столбце. Затем алгоритм исключения применяется ко всем уравнениям, кроме первого (фактически это будет i-е уравнение). На следующем шаге мы должны исключить x2 во всех уравнениях, кроме первых двух. Рассмотрим во втором столбце все элементы, начиная со второго, и выберем из них наибольший по модулю элемент. Пусть это будет элемент ak2. Поменя- ем местами второе и k-е уравнение, так что элемент ak2 станет теперь ведущим элементом, который расположен во второй строке и во вто- ром столбце. Исключим x2 во всех уравнениях, кроме первых двух. Аналогичный процесс повторяется до тех пор, пока матрица системы не примет треугольный вид. 1.2. Метод простых итераций В отличие от прямых итерационные методы решения СЛАУ дают возможность получить решение лишь приближенно на основе группы определенных, неоднократно повторяющихся операций. Однократное выполнение такой группы операций называется итерацией. Получение все более точного решения на основе многократного повторения ите- раций возможно в том случае, когда итерационный метод сходится к искомому решению. В методе итераций решение системы (1.1) строит-  16 ся в виде последовательности векторов х(0) , х(1) , х( 2) ,... , х( k ) ,... , где x  ( x1 , x2 ,..., xn ) , а верхний индекс обозначает номер итерации. При ре- шении системы (1.1) методом итераций необходимо выполнить следу- ющее: – задать начальное приближение х(0)  ( х1(0) , х2(0) ,... , хn(0) ) ; – привести исходную систему к виду, удобному для проведения ите- раций; – сформулировать последовательность действий, составляющих одну итерацию; – определить сходимость итерационного процесса к точному решению; – указать условия окончания итерационного процесса. Первый из перечисленных вопросов для систем линейных уравне- ний оказывается не очень существенным, так как, если итерационный метод сходится, то он сходится при любом начальном приближении. Поэтому обычно в качестве начального приближения берется какой-  либо известный вектор, например, вектор правых частей b или нуле- вой вектор. Конечно, если из существа решаемой задачи известно бо- лее точное приближение к решению, то надо воспользоваться им, так как в этом случае для получения решения с заданной точностью потре- буется меньше итераций. Сведение линейной системы к виду, удобному для итерации – спе- циальный способ записи линейной системы, в котором каждая искомая переменная явно выражается через все остальные путем преобразова- ния исходной системы. Для проведения итераций надо записать исход- ную линейную систему (1.1) в таком виде, чтобы на каждой итерации можно было определять уточненное значение для каждой из неизвест- ных x1, x2,... xn. Для этого надо из n уравнений (1.1) получить n явных выражений для каждой из переменных через остальные переменные  х1  1 ( х1 , х2 ,...,хn ),   х2   2 ( х1 , х2 ,...,хn ) , (1.14)  ..................................  хn   n ( х1 , х2 ,...,хn ), где 1, 1,…, n, – произвольные линейные функции. Система (1.1) может быть сведена к виду (1.14) различными способами. Стандарт- ный прием получения (1.14) заключается в том, что из первого уравне- ния выражается x1 , из второго – x 2 и т.д.  17  1  х1  а (b1  a12 x2  a13 x3 ...  a1n xn ),  11  1 (1.15)  х2  (b2  a21x2  a23 x3 ...  a2 n xn ),  а 22 .................................................................   1  хn  а (bn  an1 x1  an 2 x2 ...  an,n1 xn1 ).  nn При записи (1.15) использовано предположение, что все aii  0. Если некоторые диагональные элементы будут нулевыми, то для невы- рожденной матрицы всегда можно получить вид (1.15) простой пере- становкой строк. В общем виде метод простой итерации для системы вида (1.15) можно записать следующим образом:  (k ) 1 ( k 1) ( k 1) ( k 1)  х1  a (b1  a12 x2  a13 x3 ...  a1n xn ),  11  (k ) 1 (1.16)  х2  (b2  a21x1( k 1)  a23 x3( k 1) ...  a2 n xn( k 1) ),  a 22 .....................................................................................   (k ) 1 ( k 1) ( k 1) ( k 1)  хn  a (bn  an1 x1  an 2 x23 ...  an ,n1 xn1 ).  nn Одним из возможных условий окончания итерационного процесса является малое изменение приближенных решений от итерации к ите- рации. Это может быть, например, следующее условие: х ( k )  x ( k 1)   (i  1, n) , i i (1.17) где  – задаваемое малое число, характеризующее точность расчета. Другое возможное условие – это насколько хорошо приближенное решение x(k) удовлетворяет исходной системе (1.1) a11x1( k )  a12 x2( k ) ...  a1n xn( k )  b1  r1( k ) ,  a21x1  a22 x2 ...  a2 n xn  b2  r2 , (1.18) (k ) (k ) (k ) (k )  .................................................................. a x ( k )  a x ( k ) ...  a x ( k )  b  r ( k ).  n1 1 n2 2 nn n n n Правые части выражения (1.18) называются невязками уравнений системы. Если x (k ) является точным решением, то невязки будут рав- ны нулю. Малость невязок можно считать признаком того, что при- ближенное решение близко к точному. Аналитически данный факт означает выполнение следующего условия:  18 ri( k )  ri( k 1)   1 (i  1, n). (1.19) Для формулировки достаточного условия сходимости рассмотрен- ного метода введем понятие матрицы с диагональным преобладанием. Матрица имеет диагональное преобладание, если имеют место следу- ющие соотношения:  а11  а12  а13 ...  а1n ,   а22  а21  а23 ...  а2 n , (1.20)  ............................................  а  а  а ...  а  nn n1 n2 n ,n1. Пример 1.2. Рассмотрим две матрицы A и B: 5 2 1  3 1 1          А  3 4 1  В    4 5 2 .           8 6 15  0 1 8      Из них матрица А имеет диагональное преобладание, так как 5  2   1 , 4  3  1 , 15   8  6 , а матрица В – нет, так как во вто- рой строке 5   4  2. Условия сходимости итерационного метода – математические соотношения, при выполнении которых последовательность значений, получаемых в итерационном процессе, сходится к точному решению (в данном случае к решению системы линейных алгебраических уравне- ний). Если матрица A исходной линейной системы (1.1) имеет диаго- нальное преобладание, причем хотя бы одно из неравенств (1.20) явля- ется строгим, то метод простой итерации сходится. Пример 1.3. Требуется решить систему уравнений методом итера- ций 2 х1  х2  3х3  0,  (1.21) 4 х1  х2  х3  2, 2 х  5 х  3х  10.  1 2 3 Если записать (1.21) непосредственно в виде (1.15), выражая хi из i-го уравнения соответственно  х1  1 2  (0  x2  3x3 ) ,  (1.22)  х2  2  4 x1  x3 ,  х  1 3  (10  2 x  5 x ).  3 1 2 то метод итераций может расходиться, так как у исходной матрицы системы (1.21)  19 2  3   1   4  1  1     2 5 3    не выполняются условия диагонального преобладания 2  1   3 ,  1  4   1 , 3  2  5. Однако можно заметить, что в каждом уравне- нии (1.21) есть такая переменная хi, коэффициент при которой не меньше, чем сумма модулей коэффициентов при других неизвестных в этом уравнении. В первом уравнении – это переменная х3, во втором – х1, в третьем – х2. Если теперь поменять местами уравнения (1.21), сделав второе уравнение первым, третье уравнение – вторым, а первое – третьим 4 х1  х2  х3  2,  (1.23) 2 х1  5 х2  3х3  10, 2 х  х  3х  0.  1 2 3 то матрица этой системы 4 1  1 (1.24)    2 5 3      2 1 3    будет иметь диагональное преобладание 4   1   1 , 5  2  3 ,  3  1  2. Записав теперь систему (1.23) в виде удобном для проведе- ния итераций  х1  1 4  (2  x2  x3 ),  (1.25)  х2  1 5 (10  2 x1  3x3 ),  х  1 3  (0  2 x  x ).  3 1 2 получим сходящийся итерационный процесс, численные результаты которого приведены в табл. 1.2. Таблица 1.2 Численные результаты итерационного процесса k 0 1 2 3 4 5 (k ) 0 0,5 1,183333 1,019444 0,978426 1,001651 х1 (k ) 0 1,8 0,966667 0,925556 1,015741 1,004821 х2 (k ) 0 0,933333 1,11111 0,988148 0,990864 1,002708 х3  20 1.3. Сравнительная оценка прямых и итерационных методов Выбор численного метода решения систем линейных уравнений зависит от вида решаемой задачи и типа имеющейся вычислительной техники. Для заполненных матриц сравнительно небольшой размерно- сти (n < 1000) одним из наиболее эффективных является, по-видимому, метод Гаусса и его различные модификации, так как он прост в реали- зации и наиболее универсален. Однако с ростом n начинают прояв- ляться недостатки прямых методов, которые заключаются в следую- щем: – они требуют, как правило, хранения всей матрицы, что является очень серьезным ограничением на объем требуемой оперативной па- мяти (можно, конечно, хранить всю матрицу на жестком диске, но то- гда требуется организовать постоянный обмен массивов информации между оперативной памятью и диском, что значительно увеличивает общее время решения задачи); – происходит постоянное накапливание погрешностей округления, так как на каждом этапе вычислений используются результаты предыду- щих операций. В случае, если матрица является разреженной, но ненулевые эле- менты располагаются в матрице упорядоченным образом (например, матрица является ленточной или ненулевые элементы имеют вид кле- ток, расположенных вдоль главной диагонали), то метод Гаусса также является достаточно эффективным, так как можно организовать про- грамму реализации метода так, чтобы нулевые элементы не хранить в памяти ЭВМ. Для разреженной матрицы произвольного вида метод Гаусса ста- новится неудобным в силу того, что в процессе преобразований вместо некоторых нулевых элементов появляются ненулевые элементы, что приводит к значительному увеличению количества ненулевых элемен- тов. В этом случае итерационные методы могут быть гораздо более эффективными, так как для них, как правило, не требуется хранить нулевые элементы и совершать с ними какие-либо операции. Более того, во многих случаях, встречающихся в практических задачах, мож- но вообще не хранить элементы матрицы в памяти ЭВМ, а вычислять их по мере необходимости. Другое достоинство итерационных методов заключается в том, что в них погрешность не накапливается в процессе вычислений, так как точность вычисления каждой итерации определя- ется лишь результатами предыдущей итерации и почти не зависит от результатов других итерационных шагов. Это свойство самокорректи- ровки метода, когда допущенная на какой-либо итерации небольшая  21 ошибка не влияет на конечный результат, особенно полезно при реше- нии систем линейных уравнений высокого порядка. Недостатком итерационных методов во многих случаях является их довольно медленная сходимость, особенно при больших n. Поэтому в вычислительной математике большое внимание уделяется вопросам ускорения сходимости итерационных методов. Во многих случаях оказывается удобным комбинировать примене- ние прямых и итерационных методов. Одним из таких подходов явля- ется итерационное уточнение решения системы линейных уравнений, полученного каким-либо прямым методом. Изложим коротко суть это- го метода. Пусть имеется система линейных уравнений, записанная в матрич-   ном виде Ax  b и х* − приближенное решение, полученное каким-либо методом. Подставим х* в уравнение и вычислим невязку уравнения r*  b  Ax*. Для точного решения x невязка равна нулю r  b  Ax. Для приближенного решения невязка будет достаточно малым числом. По- строим погрешность решения 1  x*  x и получим уравнение для определения погрешности A1  A( x*  x)  Ax*  Ax  b  (b  r )  r. Приведенные выкладки показывают, что погрешность удовлетво- ряет исходному уравнению, в котором в правой части вместо вектора b стоит вектор невязки r , вычисленный подстановкой приближенного решения x в исходное уравнение. Определим из уравнения А1  r (1.26) величину  1. Если бы величина  1 определялась из (1.26) точно (без погрешностей округления), то мы определили бы точное решение х*  х  1. Однако из-за наличия погрешностей знание  1 позволяет нам определить не точное решение, а уточненное значение решения х1  х  1. (1.27) Во многих случаях вычислительный алгоритм, основанный на со- отношениях (1.26) и (1.27), позволяет существенно уточнить первона- чальное приближенное решение x при условии, что невязка r вычис- ляется достаточно точно. Вообще говоря, процесс итерационного уточнения можно продолжать, определяя величины x2, x3 ,... , однако это требует все более точного вычисления невязок уравнения. Удоб- ство использования метода итерационного уточнения заключается в том, что для определения  1 приходится решать систему линейных уравнений с той же исходной матрицей A. Поэтому, если для решения   исходной системы Ax  b используется LU-разложение, то знание  22 треугольных матриц L и U позволяет без труда решить систему (1.26) с помощью небольшого количества операций (n2 операций вместо при- мерно n3 операций в методе Гаусса). 2. НАХОЖДЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИЙ И СОБСТВЕННЫХ ВЕКТОРОВ МАТРИЦ Нахождение собственных значений и собственных векторов матриц требуется во многих физических и технических задачах при исследо- вании устойчивости различных процессов, при определении устойчи- вости и колебаний различных инженерных сооружений. В вычисли- тельной математике исследования по устойчивости алгоритмов и ана- лизу сходимости различных вычислительных методов также требуют решения задач на собственные значения. Собственные значения для матрицы n-го порядка – это n действи- тельных или комплексных чисел, удовлетворяющих (2.1) и определяе- мых как корни характеристического многочлена данной матрицы и характеризующих свойства матрицы при различных вычислительных преобразованиях. Матрица A n-го порядка имеет n собственных значе- ний, которые мы обозначим через 1 , 2 ,... , n , и соответствующие им    ненулевые собственные векторы u1 , u2 ,... , un. В общем случае собствен- ные значения могут быть как действительными, так и комплексными числами. Если среди собственных значений встречаются два или не- сколько одинаковых, то они называются кратными. В этом случае соб- ственные векторы, соответствующие кратным собственным значениям, объединяются в один вектор. По определению собственные вектора и собственные значения удовлетворяют соотношению   Aui  i ui (i  1,2,..., n) , (2.1)  где ui – ненулевые векторы. Смысл этого выражения заключается в том, что собственные векторы – это такие ненулевые векторы, которые матрица A преобразует в коллинеарные им векторы, при этом соответ- ствующие им собственные значения являются коэффициентами растя- жения (сжатия) длины исходного собственного вектора.  Обозначив координаты собственного вектора ui  {ui1 , ui 2 ,... , uin } и на основании соотношения (2.1) можно утверждать, что нахождение собственного вектора ui сводится к решению системы линейных уравнений для нахождения его координат  23 a11ui1  a12ui 2 ...  a1nuin  i ui1 ,  a21ui1  a22ui 2 ...  a2nuin  i ui 2 ,. (2.2)  ...................................................  an1ui1  an 2ui 2 ...  annuin  i uin. Система (2.2) имеет ненулевые решения лишь тогда, когда опреде- литель ее матрицы системы равен нулю и, в общем случае, можно записать условия существования ненулевого решения для произволь- ного собственного вектора в виде a11   a12... a1n a21 a11  ... a2 n (2.3)  0............. an1 an 2... ann   Раскрывая определитель (2.3), получаем характеристический мно- гочлен с0 n  с1n1  с2 n2 ...  сn1  сn  0. (2.4) Корни этого многочлена являются собственными значениями 1 , 2 ,... , n матрицы A. Для конкретного собственного значения i си- стема уравнений (2.2) позволяет найти координаты собственного век-  тора ui  {ui1 , ui 2 ,... , uin }. Отметим некоторые свойства собственных значений: 1. Если собственные значения невырожденной матрицы действи- тельны и различны, то соответствующие им собственные векторы ор- тогональны и образуют базис рассматриваемого n-мерного простран- ства. 2. Все собственные значения симметричной матрицы являются действительными числами. 3. Две матрицы A и C называются подобными, если существует та- кая невырожденная матрица P, что C = P–1A P. (2.5) Собственные значения подобных матриц совпадают. 4. Собственные значения  матрицы А и обратной ей матрицы А 1.связаны следующим соотношением: 1  , (2.6)  где  – собственное значение обратной матрицы. Нахождение всех собственных значений и собственных векторов матрицы называется полной проблемой собственных значений. Нахождение лишь некоторых из собственных значений и собственных  24 векторов называется частичной проблемой собственных значений. Во многих практических задачах приходится решать обобщенную проблему собственных значений. Она определяется следующим об- разом. Пусть заданы невырожденные матрицы A и B. Требуется найти  ненулевой вектор u и число  , удовлетворяющие равенству   Au  Bu. (2.7)  Тогда u и  соответственно называются обобщенным собствен- ным вектором и обобщенным собственным значением. Легко показать, что обобщенную задачу можно свести к обычной задаче на собственные значения и собственные векторы. Действитель- но, умножив обе части векторного равенства (4.15) слева на обратную матрицу B 1 , получим     В 1 Au  В 1 Bu  В 1 Аu  u. Видно, что обобщенная задача сводится к обычной задаче для матрицы В 1 A. Перейдем к вопросам численного нахождения собственных зна- чений и собственных векторов. Эта задача является одной из наиболее сложных вычислительных задач алгебры. На первый взгляд может по- казаться, что определение собственных значений сводится к нахожде- нию корней характеристического многочлена (2.4), однако сложность здесь заключается в том, что для произвольной матрицы (особенно в случае большой размерности) весьма затруднительно вычислить сами коэффициенты сi характеристического многочлена. Дополнительные трудности возникают из-за того, что среди собственных значений ча- сто встречаются кратные. Поэтому большинство численных методов основано не на получении характеристического многочлена (2.4) для исходной матрицы, а на различных преобразованиях этой матрицы с помощью соотношений типа (2.5) и ее сведению к более простому ви- ду, откуда собственные значения определяются гораздо проще. На практике для симметричных матриц чаще всего используются различные варианты метода вращений. Не приводя подробное описа- ние метода, которое является достаточно громоздким, опишем лишь его суть. К исходной матрице A многократно применяется преобразо- вание (2.5), в результате чего матрица сводится к гораздо более про- стому виду: к диагональной матрице или к трехдиагональной матрице. Для диагональной матрицы собственные значения могут быть вычис- лены сразу 1  а11  , 2  а22 ,..., n  аnn , а для трехдиагональной матрицы они вычисляются с помощью простых рекуррентных формул. Матрица P, с помощью которой производится преобразование (2.5), подбирается таким образом, чтобы на каждом шаге преобразова-  25 ния в преобразованной матрице P–1A P появлялись два нулевых эле- мента, симметричных относительно главной диагонали. Такой матри- цей P является матрица вращения, имеющая вид 1 0............... 0 0 1..................  (2.8) ........................  ......... cos ...  sin ...... P ........................   .........  sin ... cos ...... ........................    0 0............... 1  Эта матрица отличается от единичной матрицы E только измене- нием четырех членов: два из них стоят на главной диагонали pii = pjj = =cos , а два других симметричны относительно главной диагонали pii = pjj =sin . Геометрический смысл умножения на матрицу (2.8) заключается в том, что при таком преобразовании осуществляется плоское вращение в n-мерном пространстве в плоскости, проходящей через оси i и j на угол  (остальные оси координат при таком вращении остаются непо- движными). Угол поворота  выбирается таким образом, чтобы два элемента в исходной матрице, симметричные относительно главной диагонали, обратились в нуль. Существуют различные варианты метода вращений, которые отли- чаются тем, к какому виду сводится исходная матрица A, и способом организации вычислений. 2.1. Определение наибольшего и наименьшего собственных значений итерационным методом В практических задачах часто нужны не все собственные значения, а лишь некоторые из них. Чаще всего (например, в вопросах устойчи- вости) требуется определить лишь минимальное (или максимальное) по модулю собственное значение матрицы. В этом случае нецелесооб- разно решать полную проблему собственных значений. Для решения частичной проблемы собственных значений, состоя- щей в определении одного или нескольких собственных значений и соответствующих им собственных векторов, проще всего использовать итерационные методы. Строится итерационный процесс, который сходится к одному соб- ственному вектору с соответствующим вычислением одного собствен- ного значения. Опишем такой процесс для нахождения максимального  26 по модулю собственного значения и соответствующего собственного вектора. Рассмотрим случай, когда матрица A n-го порядка имеет все действительные собственные значения 1 , 2 ,... , n и соответствующие им собственные векторы u1 , u2 ,... , un. Пусть также все собственные значения расположены в порядке убывания, причем значение 1 по модулю строго больше, чем остальные собственные значения 1  2  2 ...  n. (2.9) Собственные значения матрицы обладают следующим свойством: если все они действительны и различны, то соответствующие им век- торы образуют базис в n-мерном пространстве. Это значит, что лю- бой n-мерный вектор можно представить в виде линейной комбинации     векторов u1 , u2 ,... , un. Возьмем произвольный вектор y 0. Тогда     n  y0  c1u1  c2 u 2 ...  cn u n   ci ui. (2.10) i 1   Умножим теперь матрицу A на вектор y 0. В результате вектор y 0  преобразуется в другой вектор y1   y1  Ay0. (2.11)  Запишем в развернутом виде вектор y1 , используя (2.10), (2.11) и определение собственных векторов и собственных значений   Aui  ui (i  1,2,... , n). (2.12) Тогда (2.11) запишется в виде   n  n  n  y1  Ay0  A( ci ui )   ci ( Aui )   ci i ui. (2.13) i 1 i 1 i 1 Формула (2.13) показывает, что действуя матрицей A на произ-  вольный вектор y 0 , мы получаем вектор, коэффициенты разложения которого по собственным векторам умножаются на соответствующее данному собственному вектору собственное значение.  Умножим теперь матрицу A на полученный вектор y1 и получим  новый вектор y 2 :    n  n  n  y2  Аy1  A2 y1  A( ci i ui )   ci i ( Aui )   ci 2i ui. i 1 i 1 i 1 Повторяя многократное умножение матрицы A на вновь получае- мые векторы, получаем    n  y k  Аy k 1  Ak y0   ci ki ui. (2.14) i 1  27 Теперь воспользуемся тем, что 1 максимальное по модулю соб- ственное значение, и преобразуем (2.14)         y k  1k [c1u1  c2 ( 2 ) k u 2  c3 ( 3 ) k u3 ...  cn ( n ) k u n ]. 1 1 1 Так как ввиду (2.9) 2    1, 3  1,..., n  1, 1 1 1 то с ростом степени k эти величины убывают и при k  все (i 1 ) k 0. Это означает, что при достаточно больших k коэффици-    енты, на которые умножаются векторы u2 , u3 ,... , un , будут малы и, сле-  довательно, вектор y k будет почти коллинеарен первому собственно- му вектору   yk  ki c1u1. (2.15) Приведенные соотношения (2.14), (2.15) являются основой степен- ного метода нахождения максимального собственного значения. Сте- пенной метод – итерационный метод нахождения наибольшего по мо- дулю собственного значения матрицы и соответствующего ему соб- ственного вектора. Метод называется степенным, так как в нем вычис- ляются различные степени матрицы. Однако алгоритм его реализации на практике несколько отличается от того, который был описан выше. Это вызвано тем, что при достаточно большом 1 длина получаемых по формуле (2.14) векторов сильно растет. Поэтому удобно реализо- вать степенной метод таким образом, чтобы на каждом шаге итерации вида (2.14) проводить нормировку получаемого вектора, делая его длину, равной единице (здесь мы пользуемся тем, что собственный вектор определяется с точностью до произвольного сомножителя и поэтому, разделив каждую компоненту вектора на его длину, мы полу- чим вектор единичной длины). Один из удобных способов реализации степенного метода заключается в следующем.  Задается начальный вектор y 0 с координатами { y01, y02 ,... , y0n }. Вы-  числение следующего приближения – вектора y1 – состоит из трех этапов.  1) вычисляется длина вектора y 0 n d0  ( y i 1 0i )2 ; (2.16)  2) производится нормировка вектора y 0 путем деления каждой его компоненты на величину d0  28   W0  1 y0. (2.17) d0  В результате вектор W0 будет иметь длину, равную единице;  3) вычисляется вектор y1 по формуле   y1  AW0. (2.18) На этом завершается одна итерация степенного метода.   В общем случае переход от y k к вектору y k 1 совершается по формулам n dk  ( yi 1 ki )2 ; (2.19)   Wk  1 yk ; (2.20) dk   yk 1  AWk. (2.21) Итерации проводятся до тех пор, пока не выполнится условие d k 1  d k  .   Тогда можно принять, что u1  yk , 1  d k. Это означает, что после- довательность d k сходится к максимальному собственному значению   1 , а последовательность векторов y k – к собственному вектору u1 , соответствующему значению 1. Пример 2.1. Найти максимальное собственное значение и соответ- ствующий собственный вектор для матрицы A 10 1 . А     1 2  Решение. Зада?