Размещено на http://www.allbest.ru/

БЕЛОРУССКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

ФАКУЛЬТЕТ РАДИОФИЗИКИ И ЭЛЕКТРОНИКИ

С. Г. Мулярчик

Численные методы

Конспект лекций

Минск - 2001

Содержание

Введение

Лекция 1. Введение в теорию численных методов

1.1 Основные требования, предъявляемые к вычислительным алгоритмам

Лекция 2. Системы линейных алгебраических уравнений

Лекция 3. Устойчивость и точность прямых методов

3.1 Устойчивость алгебраических методов

3.2 Точность алгебраических методов

3.3 Устойчивость метода Гаусса

3.4 Устойчивость и точность метода LU-факторизации

Лекция 4. Итерационные методы

4.1 Схема итерационных методов

4.2 Метод простой итерации (Метод Якоби)

4.3 Метод последовательной релаксации

Лекция 5. Модификации метода сопряженных градиентов

5.1 Основные свойства метода сопряженных градиентов

5.2 Экономичная форма метода сопряженных градиентов

5.3 Обобщенный метод сопряженных градиентов

Лекция 6. Численное решение уравнений. Метод итераций

6.1 Теорема о сходимости и точности метода итераций

Лекция 7. Метод ньютона для систем и его модификации

7.1 Метод Ньютона для решения систем нелинейных уравнений

7.2 Модификации метода Ньютона. Краткий обзор

7.3 Дискретный метод Ньютона

Лекция 8. Полиномиальная интерполяция

8.1 Интерполяционная формула Лагранжа

8.2 Интерполяционная формула Ньютона

Лекция 9. Сплайн-интерполяция

9.1 Ограничения многочленной интерполяции

Лекция 10. Методы численного интегрирования

10.1 Формула Симпсона для вычисления двойных интегралов

10.2 Вычисление многомерных интегралов методом Монте-Карло

Лекция 11. Численное решение систем обыкновенных дифференциальных уравнений

11.1 Общая характеристика численных методов

Лекция 12. Одношаговые методы

12.1 Линейные многошаговые методы

Лекция 13. Многошаговые методы

13.1 Явные многошаговые методы Адамса

13.2 Неявные многошаговые методы Адамса

Лекция 14. Устойчивость многошаговых методов

14.1 Построение области устойчивости

14.2 Начальные вычисления в многошаговых методах

Лекция 15. Изменение шага в многошаговых методах

Литература

Введение

В курсе лекций Численные методы рассматриваются компьютерные методы решения типовых задач математики, которые наиболее часто встречаются на практике. К ним относятся: приближение функций, вычисление определенных интегралов, решение систем линейных алгебраических уравнений, решение систем нелинейных алгебраических уравнений, решение систем обыкновен-ных дифференциальных уравнений.

Данная дисциплина читается для студентов специальностей Н.02.02.00 - Радиофизика, Н.02.03.00 - Физическая электроника. В учебных планах этих специальностей имеется курс Методы вычислительного эксперимента, в котором изучаются численные методы решения задач математической физики, а также численные методы оптимизации. Поэтому эти задачи не вошли в настоящий курс.

Численные методы по своей сути являются приближенными. Один из основных вопросов, на который следует ответить при разработке численного метода - какова точность получаемого резуль-тата Другой, не менее важный, который порождает опять-таки работа в условиях ошибок - какова устойчивость метода к возникающим на операциях ошибкам Важные характеристики метода - его эффек-тивность по времени решения задачи и экономичность по требуемым для этого объемам памяти. Ответы на эти вопросы по каждому классу рассматриваемых математических задач можно найти в лекциях.В вычислительной математике для перечисленных выше математических задач разработано большое количество численных методов. Отбирая материал для той или иной лекции, автор исполь-зовал свой опыт решения задач компьютерного моделирования. В лек-ции включены только базовые алгоритмы, лежащие в основе многих численных методов; алгоритмы, популярные у начинающих вычис-лителей; наконец, алгоритмы, обладающие высокой эффективностью.

Лекция 1. Введение в теорию численных методов

ЦЕЛЬ ЛЕКЦИИ: определить предмет курса; дать общую характеристику таких свойств численных методов как эффективность, точность, устойчивость, экономичность; пояснить эти свойства на простейших вычислительных правилах.

1.1 Основные требования, предъявляемые к вычислительным алгоритмам

Вычислительный алгоритм должен быть устойчивым к ошибкам, обеспечивать требуемую точность, быть эффективным по времени выполнения, быть экономичным по требуемым объемам памяти, не допускать аварийных остановов. Поясним эти требования.

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

при заданных и d.

Предположим, что при вычислении внесена ошибка i (например, за счет операции округления), т. е. вместо значения использу-ется приближенное значение

.

В соответствии с рекурсивным правилом

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

и ошибка, допущенная на i-м шаге процесса, на следующем шаге не увеличивается, если операция сложения выполнена без новых округлений. Это означает, что алгоритм устойчив.

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

Опять y0 и q считаем заданными.

Пусть, как и в предыдущем примере,

.

Тогда

,

или

.

В этом случае

,

и при q>1 ошибка будет возрастать. Такой алгоритм является неустойчивым.

Точность. Обозначим через y точное решение задачи, yk - ее при-ближенное значение, полученное на k-м шаге численного метода. Тогда

составляет погрешность решения, а

является абсолютной погрешностью.

Вычислительный алгоритм должен давать решение с заданной точностью и, следовательно, критерием завершения процесса уточнения решения является выполнение неравенства

.

На практике в силу трудностей вычисления абсолютной погрешности вместо используют ее оценку сверху - предельную абсолютную погрешность Дy. В качестве Дy выбирают как можно меньшее значение, удовлетворяющее неравенству

,

а критерием завершения процесса в этом случае является неравенство

.

Часто используют относительную погрешность



т. е. абсолютную погрешность на единицу измерения, и предельную относительную погрешность

.

Критерий завершения процесса в этом случае имеет вид

,

где - заданная допустимая относительная ошибка.

Рассмотрим в качестве примера, как может быть построена оценка предельной абсолютной погрешности вычисления значений функции

,

если известны абсолютные погрешности ее аргументов.

Ошибка имеет вид

.

Тогда

,

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



Пусть решением задачи является вектор

,

который будем рассматривать как элемент векторного пространства . Приближенное решение

и его погрешность

также являются элементами этого пространства. Рассмотрим, как оцениваются ошибки решения в этом случае. Для этого используем количественную характеристику вектора в виде нормы.

Говорят, что в задана норма, если каждому вектору из сопоставляется вещественное число, называемое нормой вектора , для которого справедливы следующие свойства:

, причем тогда и только тогда, когда ,

для любого вектора и любого числа ,

для любых векторов и .

В вычислительных методах наиболее употребительны следующие нормы:

.

Абсолютную и относительную погрешность вектора в любой из перечисленных выше норм можно определить следующим образом:

.

Имеет смысл говорить об абсолютной и относительной погрешности (mn)- матрицы решений A. В этом случае используется такая характеристика, как норма матрицы, согласованная с нормой вектора. Для нормы матрицы A, обозначаемой , справедливы следующие свойства:

, причем тогда и только тогда, когда,

для любой матрицы A и любого числа ,

для любых (mn)-матриц и.

Каждой из векторных норм соответствует своя согласованная норма матриц. В частности, нормам соответствуют нормы , вычисляемые по правилам

где j(ATA) - собственные числа матрицы ATA.

Абсолютная и относительная погрешности матрицы решений вычисляются через соответствующие нормы:

,

причем - искомая матрица решений, - ее некоторое приближение.

Эффективность. Эффективность вычислительного алгоритма оценивают по времени реализации либо по общему числу операций, требуемых для получения результата.

Обозначим через nсл число операций сложения, через nу - число операций умножения, через nдел - число операций деления, через tу, tдел, tсл - время выполнения операций умножения, деления и сложения на ЭВМ. Объективный критерий сравнения вычислительных алгоритмов - требуемое время вычислений:

T=nуtу+nслtсл+nделtдел.

Если tу?tдел?tсл, то T?(nу+nдел+nсл)tоп, где tоп - время выполнения на ЭВМ арифметической операции, и критерием оценки эффективности алгоритма может быть

nУ=nу+nдел+nсл .

Если tу?tдел>>tсл, то критерий эффективности - число так называемых длинных операций nдл=nдел+nу.

Пусть Q - критерий эффективности метода. Значение Q зависит от требуемой точности . Из двух сравниваемых вычислительных алгоритмов решения задачи из них будет эффективнее тот, который при заданном значении е характеризуется меньшим значением Q.

В качестве примера построим два алгоритма вычисления полинома n-й степени

.

Непосредственная реализация:



Расчет полинома по рекурсивному правилу

,

вытекающему из приведенной реализации, соответствует алгоритму 1.

Алгоритм 1.

Начальная установка

Вычислить в цикле (k=1,2,…,n):

Алгоритм 1 характеризуется таким числом арифметических операций: nу=2n; nсл=n.

Преобразуем исходную запись полинома:

Теперь расчет полинома можно выполнить следующим образом:

Записанное правило реализуется алгоритмом 2.

Алгоритм 2.

Начальная установка:

Вычислить в цикле (k=n,n-1,…,1):

Значение . Затраты на выполнение алгоритма:

nу=n , nсл=n.

Очевидно, что алгоритм 2 эффективнее по времени выполнения алгоритма 1.

Экономичность. Экономичность алгоритма оценивается по требуемым для его реализации объемам памяти ЭВМ.

Аварийные остановы. Любое число ЭВМ принадлежит не всей числовой оси, а некоторому интервалу . Если какое-либо число x в процессе работы алгоритма выходит за указанный интервал, то происходит так называемый аварийный останов. Необходимо построить вычислительный алгоритм таким образом, чтобы вычисленные в процессе его работы числа .

Лекция 2. Системы линейных алгебраических уравнений

ПРЯМЫЕ МЕТОДЫ

ЦЕЛЬ ЛЕКЦИИ: Определить два класса численных методов (прямые и итерационные); показать, как строятся прямые методы Гаусса, LU-факторизации, Холесского; выполнить оценку их эффективности.

Постановка задачи.

Основная задача вычислительной алгебры - решение систем линейных алгебраических уравнений (СЛАУ)

В дальнейшем будем использовать запись этой системы в компактной форме:

( запись означает, что индекс i изменяется от 1 до n с шагом 1), или в векторном виде

,

где

Предполагается, что матрица неособенная, т. е. , и решение единственно.

Прямые и итерационные методы.

Численные методы решения СЛАУ делятся на две большие группы: прямые и итерационные.

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

,

где k - номер итерации. В действительности итерационный процесс прекращается, как только становится достаточно близким к .

Имеется промежуточный класс методов, в которых решение ищется итерационно, однако для них заранее известно, какое число итераций необходимо выполнить, чтобы в отсутствии ошибок округления получить точное решение. На практике при вычислении приближенного решения число итераций в наиболее эффективных методах оказывается значительно меньшим, чем этого требует теория точного решения.

Какой класс методов лучше? Однозначно на этот вопрос ответить нельзя. Итерационные методы привлекательнее с точки зрения объема вычислений и требуемой памяти, когда решаются системы с матрицами высокой размерности. При небольших порядках системы используют прямые методы либо прямые методы в сочетании с итерационными методами.

Метод Гаусса.

В методе Гаусса линейная система

решается в два этапа. На первом этапе система преобразуется к виду ,



Рис. 2.1. Структура системы и портрет ее ненулевых элементов до (а) и после (б) прямого хода Гаусса

где - верхняя треугольная матрица с единичной диагональю (это так

называемый прямой ход Гаусса). На втором этапе (обратный ход Гаусса) решается система . Рассмотрим эти этапы подробнее.

Прямой ход. Прямой ход Гаусса состоит из n шагов.

Первый шаг. Полагаем, что и разделим на него первое уравнение. Перепишем систему с учетом этого преобразования:

Умножим первое уравнение на и вычтем его из i-го уравнения преобразованной системы:



Обозначим

.

Получим

Второй шаг. На втором шаге из системы

исключается аналогичным образом:

K-й шаг. Запишем общий вид преобразованной системы после k-го шага прямого хода Гаусса:

Здесь

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

Оценим количество длинных операций (умножений и делений) на первом шаге прямого хода Гаусса. Преобразование первого уравнения требует n таких операций. Преобразование остальных n-1 уравнений - n(n-1) операций умножения и деления. Таким образом, первый шаг выполняется за длинных операций. Рассуждая по аналогии, нетрудно найти затраты на остальных n-1 шагах. Суммарные затраты прямого хода Гаусса определяются в итоге рядом

.

Последняя оценка имеет место для n>>1.

Обратный ход. Запишем систему, решаемую на обратном ходе, в координатном виде

Ее решение:



Запись означает, что индекс k изменяется от значения n-1 до 1 с шагом 1.

Требуемое число длинных операций на обратном ходе

Приближенная оценка справедлива для n>>1.

Общие затраты метода Гаусса:

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

Метод LU-факторизации.

В методе LU-факторизации (эту схему называют компактной схемой Гаусса) при решении системы выполняется следующая последовательность действий.

Матрица представляется в виде произведения

,

где L нижняя треугольная матрица, U верхняя треугольная матрица. Такое разложение единственно при условии предварительного выбора диагональных элементов одной из матриц. В этом случае число элементов в матрице A совпадает с суммарным числом неизвестных элементов матриц L и U. Если диагональ L принимается единичной, то такое разложение называют разложением Дулиттла если единична диагональ U - разложением.

Система заменяется системой

,

легко решаемой за два шага:

Шаг 1. . Принимая во внимание треугольный вид матрицы L, нетрудно получить, что в алгоритме Краута

Шаг 2. . Решение этой системы в алгоритме Краута:

.

Суммарные затраты реализации обоих шагов при n>>1 составляют длинных операций.

Получим соотношения для расчета элементов матриц L и U в алгоритме Краута. Для этого перемножим матрицы L и U и приравняем результат к A. По правилу перемножения матриц

Учтем, что

Рассмотрим элемент (рис. 2.4), расположенный на центральной диагонали либо в нижней треугольной части матрицы A. Для такого элемента i ? j. Из рисунка следует, что

Иллюстрация вычисления элемента матрицы, расположенного ниже главной диагонали

,

так как i ? j и . Отсюда

Рассмотрим элемент (рис. 2.5), находящийся выше главной

Иллюстрация вычисления элемента матрицы, расположенного выше главной диагонали диагонали матрицы A (для него j>i). В этом случае.

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

Получили в итоге соотношения, которые позволяют вычислять элементы матриц L и U. Последовательность вычислений: сначала вычисляется столбец матрицы L, далее строка матрицы U, затем опять столбец матрицы L, далее строка матрицы U и т. д. (см. рис. 2.6, который иллюстрирует последовательность вычислений и схему хранения матриц L и U).

Вычисление столбца матрицы L и строки матрицы U назовем шагом LU-разложения. Приведем в качестве примера схему хранения элементов матриц A,L,U после второго шага LU-разложения (рис. 2.7).

Число длинных арифметических операций на этапе LU-разложения при n>>1 составляет величину , на шаге решения ли-

нейных систем с треугольными матрицами - . Суммарное число длинных операций приближенно равно (как и в методе Гаусса),

Рис. 2.6. Исходная матрица A (а), схема хранения L и U матриц (б), последовательность вычисления элементов в принятой схеме хранения (в)

Рис. 2.7 Схема хранения элементов 44 матриц A, L, U после второго шага LU-факторизации

т. е. основные затраты приходятся на LU-факторизацию матрицы A. Эта особенность делает особо привлекательным метод LU-факторизации при решении СЛАУ с одной и той же матрицей A, но разными правыми частями:

В этом случае факторизация матрицы выполняется однократно, требуя длинных операций, а решение каждой системы с соответствующей правой частью реализуется за таких операций.

Метод Холесского.

Исключительно эффективную реализацию метода LU-факторизации можно получить, если ограничиться классом линейных систем с симметрической положительно определенной матрицей A, т. е. Такую реализацию называют методом Холесского, либо методом квадратного корня.

Будем полагать, что решаемая система

имеет симметрическую положительно определенную матрицу A. В этом случае матрица A представляется в виде

Здесь - нижняя треугольная матрица. Такое разложение существует и единственно для положительно определенных симметрических матриц.

Система преобразуется к виду .

Вектор ищется путем последовательного решения двух систем с треугольными матрицами: ; .Для получения расчетных соотношений элементов матрицы рассмотрим произвольный элемент матрицы A:

Суммирование здесь выполняется только до j, т. к. j?i. Выделим член при значении k=j:

.

Теперь

Лекция 3. Устойчивость и точность прямых методов

ЦЕЛЬ ЛЕКЦИИ: Выполнить оценку устойчивости и точности прямых методов; показать, как перестановкой строк и столбцов обеспечивается устойчивость и точность прямых методов, каким образом осуществляется выбор матриц перестановки строк и столбцов в случае систем с разреженной матрицей.

3.1 Устойчивость алгебраических методов

Прямые методы в отсутствии ошибок округления позволяют получить точное решение системы .

Современные вычислительные машины оперируют с конечными десятичными дробями, представленными в форме с плавающей точкой. В этом случае уже на этапе запоминания элементов матрицы A и вектора в памяти ЭВМ вносится ошибка округления и реально решается возмущенная система

.

Здесь и - возмущения матрицы системы и вектора правой части.

Для элементов матрицы и компонент вектора справедливы оценки

,

где- элемент матрицы A, - компонент вектора , е - число, характеризующее относительную погрешность машинной ариф-метики. Например, для двоичных ЭВМ, использующих арифметику с плавающей точкой и t - разрядную мантиссу, .

Перейдем к более общей числовой оценке возмущений - норме. Из записанных выше неравенств следует, что

,

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

Поясним теперь суть одного из наиболее разработанных подходов к анализу устойчивости алгебраических методов.

Пусть - приближенное решение СЛАУ, полученное применением некоторого прямого метода. Очевидно, что вследствие ошибок округления при реализации на ЭВМ прямого метода это решение не будет точно удовлетворять системе . Пусть, однако, удается показать, что является точным решением системы

.

Если для матрицы F и вектора , называемых эквивалентными возмущениями метода, можно получить оценки вида

,

где f(n), h(n) некоторые степенные функции типа с небольшим показателем k, то метод считается устойчивым по Уилкинсону. Такой вид функций f(n), h(n) объясняется тем, что в процессе реализации прямых методов на ЭВМ неизбежно некоторое накопление ошибок округления, пропорциональное числу арифметических операций, за которое прямой метод приводит к решению.

3.2 Точность алгебраических методов

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

,

где - точное решение, в любой согласованной норме подчиняется неравенству

.

Обозначим

.

Величины , являются оценками относительных возмущений матрицы A и вектора b. Число ч есть не что иное, как число обусловленности матрицы A.

Естественно допустить, что

.

Такое неравенство означает, что норма матрицы возмущений существенно меньше нормы исходной матрицы, т. е. предполагается, что прямой метод решения линейной системы устойчив. В этом случае

.

Относительная ошибка решения системы прямым методом может достигать величины, определяемой произведением числа обусловленности и суммы относительных возмущений матрицы системы и вектора правой части. Если число обусловленности велико, то даже при малых эквивалентных возмущениях матрицы и вектора правой части погрешность решения линейной системы оказывается значительной.

Задачи, для которых ч велико, называют плохо обусловленными. Число обусловленности зависит от выбора матричной нормы. Например, для симметричной положительно определенной матрицы A в качестве можно взять ее максимальное собственное число. Тогда число обусловленности

,

где , - максимальное и минимальное собственные числа матрицы.

3.3 Устойчивость метода Гаусса

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

.

Запишем оценку нормы матрицы возмущения:

.

Вид этой оценки удовлетворял бы критерию устойчивости Уилкинсона, если бы множитель g(A) имел небольшое значение. Поясним смысл множителя g(A).

Пусть обозначает матрицу, полученную из A после k шагов исключения. Обозначим

.

Тогда

.

Следовательно, g(A) показывает, во сколько раз могут возрасти элементы матрицы A в ходе исключения переменных по сравнению с их исходным уровнем. По этой причине g(A) называют коэффициентом роста матрицы A .

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

.

Для ограничения роста элементов матрицы в процессе гауссова исключения желательно, чтобы поправочные члены

в этой формуле были не слишком большими. Это достигается процедурой выбора элемента , который называют главным.

Выбор главного элемента по столбцу. В этом случае ограничение роста элементов матрицы Ak на k-м шаге гауссова исключения достигается перестановкой строк таким образом, чтобы гарантировать неравенство

.

С этой целью при исключении переменной в качестве главного элемента выбирается элемент матрицы Ak-1 по правилу

,

т. е. наибольший по модулю элемент в k-м столбце матрицы Ak-1

Строки r и k переставляются и только после этого выполняется k-й шаг исключения прямого хода Гаусса.

При столбцовой стратегии выбора главных элементов справедлива такая оценка для значения параметра ak, определяющего коэффициент роста:

.

Она допускает, что и, следовательно, коэффициент роста

.

По этой причине метод Гаусса с выбором главного элемента по столбцам является условно устойчивым. Несмотря на это, он широко используется на практике, так как g(A) редко достигает своего верхнего предела.

Выбор главного элемента по всей матрице. В этой стратегии в качестве главного элемента при исключении неизвестной xk выбирается элемент по правилу

,

т. е. наибольший по модулю элемент в квадратной подматрице матрицы Ak-1 (рис. 3.2). Строки k и r, а также столбцы k и l переставляются и далее выполняется k-й шаг исключения. Такая стратегия гарантирует выполнение неравенства

и, следовательно, ограничивает рост элементов в процессе исключения Гаусса.

Оценка коэффициента роста элементов матрицы A в этом случае имеет более благоприятный вид:

.

Точность метода Гаусса.

Привлекая оценку нормы матрицы возмущения, можно записать, что

.

Анализ неравенства позволяет определить пути повышения точности метода Гаусса: выбор главных элементов, работа с числами удвоенной длины, переобусловливание системы линейных алгебраических уравнений.

3.4 Устойчивость и точность метода LU-факторизации

Метод LU-факторизации, если не делать никаких специальных усилий, характеризуется такими же оценками нормы матрицы возмущения и относительной ошибки, как и метод Гаусса. В то же время этот метод наряду с упоминавшимися при рассмотрении метода Гаусса направлениями повышения точности решения имеет еще одну возможность: использование так называемой операции скалярного накопления. Эта операция во многих ЭВМ реализована аппаратно либо моделируется специальными подпрограммами. Операция скалярного накопления привлекается при вычислении скалярного произведения

.

Пусть основным режимом работы ЭВМ является арифметика с плавающей точкой и t-разрядной мантиссой машинного слова. Рассмотрим вычисление скалярного произведения. Произведение чисел вk и гk, выполненное арифметическим процессором, имеет 2t-разрядную мантиссу. Вместо того, чтобы округлить это произведение до t разрядов и прибавить результат к хранимой t-разрядной промежуточной сумме Sk-1, режим накопления предусматривает удержание 2t разрядов как для Sk-1, так и для вkk. Очередное суммирование выполняется с 2t-разрядными числами. Оно также сопровождается ошибкой округления, но эта ошибка затрагивает лишь последний разряд 2t-разрядного результата. Возвращение к t-разрядному представлению происходит лишь после вычисления скалярного произведения.

Анализ расчетных соотношений метода LU-факторизации свидетельствует о том, что они реализуются с использованием операции скалярного накопления.

Оценка нормы матрицы возмущений метода LU-факторизации при использовании операции скалярного накопления имеет вид

.

Устойчивость вычислительной схемы метода определяется величиной коэффициента роста элементов матрицы A. Чтобы уменьшить величину g(A), в методе LU-факторизации используется столбцовая процедура выбора главных элементов.

Погрешность решения метода

.

Устойчивость и точность метода Холесского.

Этот метод применяется для решения линейных систем с симметрической положительно определенной матрицей A. Можно показать, что в схеме Холесского коэффициент роста элементов матрицы A равен 1. Действительно, из соотношения

следует, что . Если учесть, что в симметрической положительно определенной матрице A наибольший по модулю элемент положителен и находится на главной диагонали, нетрудно прийти к сформулированному заключению: g(A)=1.

Таким образом, чтобы ошибки округления при реализации схемы Холесского были малыми, не требуется никаких перестановок строк либо столбцов. Наименьшая погрешность численного решения достигается при привлечении операции скалярного накопления. В этом случае

.

Прямые методы для систем с разреженными матрицами.

Во многих приложениях приходится решать систему

с разреженной матрицей A. Назовем -матрицу A разреженной, если число ее ненулевых элементов .

При решении таких систем естественно хранить в памяти ЭВМ только ненулевые элементы. Можно надеяться, что в случае разреженных матриц как вычислительная работа, так и требуемый объем памяти значительно сократятся. Проблема состоит в том, что на этапе исключения неизвестных в методе Гаусса либо на этапе построения L и U матриц в компактной схеме Гаусса могут возникать

Рис. 3.3. LU-факторизация диагональной матрицы со строчным и столбцовым окаймлением

новые ненулевые элементы. Обратимся к следующему простому при-меру. Хотя исходная матрица A здесь имеет ничтожное число ненулевых членов, матрицы L и U имеют такую же структуру ненулевых элементов, как и для полностью заполненной матрицы. В

то же время перестановкой строк и столбцов систему можно преобразовать к виду, который сохраняет свойство разреженности исходной матрицы в матрицах разложения. В общем случае, конечно, не всегда удается полностью сохранить разреженность. Поэтому перестановки строк и столбцов выбирают такими, чтобы заполняемость матрицы в процессе исключения переменных была минимальной.

На языке матричного описания эта задача выглядит следующим образом. Пусть P и Q - -матрицы перестановок строк и столбцов, причем . Матрицы P и Q легко получить из единичной матрицы, в которой переставлены соответствующие строки (столбцы). Например, матрица перестановки первой и третьей строк -матрицы приведена на.

Преобразуем систему линейных алге-браических уравнений следующим образом:

где - переупорядоченная форма A, , . К системе будем применять прямой метод решения (его эффективность зависит от выбора матриц P и Q) и далее найдем , так как . Задача осложняется тем, что перестановка строк и столбцов определяющим образом влияет не только на разреженность матриц разложения, но и на численную устойчивость прямого метода.

Симметричные матрицы. Если матрица A линейной системы является симметрической положительно определенной, то процесс исключения переменных по схеме Холесского является численно устойчивым при любом выборе главных элементов из числа диагональных членов. Таким образом, если положить , т. е. при перестановке k-й и j-й строк осуществлять перестановку k-го и j-го столбцов, то выбор матрицы P можно производить только исходя из требования минимизации заполнения матрицы .

Задача численного решения линейной системы сводится в этом случае к выполнению трех последовательных шагов:

Шаг 1. Формирование матрицы и вектора .

Шаг 2. Решение упорядоченной системы .

Шаг 3. Вычисление вектора .

Рассмотрим более подробно, как реализуется Шаг 1. Вообще говоря, для -матрицы A существует различных упоря-дочиваний, из которых одно или более оказывается оптимальным в смысле заполнения. Задача отыскания оптимальных упорядочиваний является чрезвычайно трудной и практически нерешенной на сегодняшний день. Существующие процедуры реализации шага 1 являются эвристическими.

Они обеспечивают понижение заполнения матрицы , однако не гарантируют достижение точного минимума. Опишем в общих чертах одну из таких процедур - алгоритм минимальной степени. Основная идея алгоритма минимальной степени заключается в локальной минимизации заполнения за счет выбора главного элемента в той строке и столбце, которые обеспечивают внесение наименьшего числа ненулевых элементов в множитель . Рассмотрим, например, k-й шаг исключения прямого хода Гаусса. В поддиагональных позициях столбцов с 1-го по (k-1)-й расположены нулевые элементы.

Для исключения переменной xk на этом шаге k-я строка, нормированная соответствующим образом, вычитается из всех тех строк, которые имеют ненулевые элементы в k-м столбце. Следовательно, чтобы минимизировать число новых ненулевых элементов, достаточно просмотреть активную подматрицу матрицы Ak-1, выбрать строку с наименьшим числом ненулевых элементов, назовем ее l-й строкой, и переставить как строки, так и столбцы с номерами l, k. После такой перестановки в матрицу Ak-1 вносятся новые ненулевые элементы, т. е. формируется портрет матрицы Ak , и процесс имитации исключения переменных продолжается. Выполнив имитацию (n-1)-го шага исключения, построим в итоге из единичной матрицы I искомую матрицу перестановок P.

Матрицы общего вида. На k-м шаге алгоритма Гаусса при решении линейных систем с разреженной матрицей A произвольной структуры необходимо выбирать в качестве главного элемента такой элемент , который прежде всего обеспечивает устойчивость вычислительного процесса и в то же время по возможности минимизирует заполняемость матрицы Ak. Это достигается введением специального барьера - числа u?1 и выделением множества элементов из активной подматрицы матрицы Ak-1, которые удовлетворяют условию

,

где - элемент множества Sk-1.

Для каждого из таких элементов вводится числовая оценка

,

где r(l,k) - число ненулевых элементов в l-й строке активной части матрицы Ak-1, c(m,k) - число ненулевых элементов в m-м столбце этой же части матрицы Ak-1. В качестве главного элемента на k-м шаге прямого хода Гаусса выбирается такой элемент

,

для которого

.

Лекция 4. Итерационные методы

ЦЕЛЬ ЛЕКЦИИ: Дать общую схему линейных итерационных методов первого порядка и построить на её основе методы простой итерации, Гаусса-Зейделя, последовательной релаксации; сравнить эти методы по вычислительным затратам на итерации.

В случае больших разреженных СЛАУ предпочтение отдаётся итерационным методам, т. к. они, во-первых, не приводят к появлению на итерациях новых ненулевых элементов, во-вторых, оказываются более эффективными по затратам машинного времени.

4.1 Схема итерационных методов

Пусть требуется решить систему

.

Предположим, что Q - неособенная -матрица, рассматриваемая как апроксимация матрицы A, для которой решение системы

не представляет особого труда. Это имеет место тогда, когда матрица Q является нижней треугольной, верхней треугольной, трёх-диагональной, произведением конечного числа таких простых матриц и т.п. Матрицу Q называют матрицей расщепления.

Перепишем исходную систему в виде

.

Построим итерационный процесс



Разрешим его относительно:

,

или

.

Введем обозначения:

; .

Тогда итерационную схему можно представить в виде

Это линейная схема первого порядка. Поскольку итерационная матрица B на итерациях является постоянной, то такую схему называют стационарной.

Достаточное условие сходимости этой итерационной схемы имеет вид .

Убедимся в этом.

Обозначим вектор погрешности на k-м шаге через , на (k+1)-ом шаге -, точное решение - . Тогда

.

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

Отсюда

.

При видно, что , т. е. метод сходится к решению. Построим теперь вычислительные схемы, выбирая конкретный способ формирования матрицы раcщепления.

4.2 Метод простой итерации (Метод Якоби)

Пусть требуется решить систему

.

Представим , где D - диагональная матрица с диагональными членами матрицы A; - часть матрицы A, лежащая ниже центральной диагонали, - часть матрицы A, лежащая выше центральной диагонали. Тогда

,

или

.

Запишем итерационный метод

Разрешим его относительно :

Эта форма удобна для реализации метода Якоби. Здесь итерационная матрица

.

Матрица расщепления . Следовательно, при построении метода Якоби в качестве матрицы ращепления используется диагональ матрицы СЛАУ.

Запишем метод Якоби в координатной форме для системы общего вида:

Нетрудно убедиться, что метод Якоби в координатной форме

есть не что иное, как разрешение каждого из уравнений системы относительно одной из компонент вектора. Из первого уравнения системы выражается и его значение принимается за значение . Из второго уравнения определяется и его значение принимается за и т. д. Переменные в правой части этих соотношений при этом полагаются равными их значениям на предыдущей итерации.

Рассмотрим пример системы линейных алгебраических уравнений с разреженной (пятидиагональной) матрицей A и оценим эффективность решения такой системы методом Якоби. Будем хранить матрицу A в виде пяти одномерных массивов g, f, a, b, c. Тогда решаемая система в терминах указанных массивов примет вид

.

В этой записи следует учесть, что отдельные коэффициенты являются нулевыми. В частности

при i=1; при ; при

;

при

Рис. 4.1. Система с пятидиагональной матрицей

Реализация метода Якоби для такой системы линейных алгебраических уравнений осуществляется просто:



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

Метод Гаусса-Зейделя.

Решаемую систему запишем в виде

.

Итерационная схема Гаусса-Зейделя также следует из этого пред-ставления системы:

,

или

.

Последняя форма такого итерационного метода удобна для реа-лизации.

Приведем метод Гаусса-Зейделя к стандартному виду:

Стандартная форма метода позволяет выписать его итерационную матрицу и провести над ней очевидные преобразования:

.

Матрица расщепления

cодержит всю нижнюю треугольную часть матрицы A.

Запишем метод Гаусса-Зейделя в координатной форме для системы общего вида:

Координатная форма метода Гаусса-Зейделя отличается от координатной формы метода Якоби лишь тем, что первая сумма в правой части итерационной формулы содержит компоненты вектора решения не на k-й, а на (k+1)-й итерации. Но именно это означает, что матрица расщепления метода Гаусса-Зейделя включает не только диагональные члены матрицы A, но и все члены, расположенные ниже главной диагонали.

Оценим, как и ранее, временные затраты метода при решении тестовой задачи СЛАУ с пятидиагональной матрицей. Реализация метода Гаусса-Зейделя для такой системы:

позволяет констатировать, что временные затраты на итерации метода



оказываются такими же, как и в методе Якоби. Ключ к пониманию более высокой эффективности метода Гаусса-Зейделя количество итераций, требующееся для получения решения с заданной точностью. Более высокая сходимость к решению в методе Гаусса-Зейделя достигается за счет выбора матрицы расщепления, лучше аппроксимирующей матрицу A.

4.3 Метод последовательной релаксации

Воспользуемся методом Гаусса-Зейделя, записанным в форме

,

или

.

Удобная для реализации итерационная схема релаксационного метода имеет вид

,

где щ - некоторый параметр.

Преобразуем ее к стандартной форме:

.

В свою очередь итерационная матрица

.

В итоге матрица расщепления

.

При указанная итерационная схема соответствует методу последовательной верхней релаксации, при - методу последовательной нижней релаксации.

Координатная форма метода последовательной релаксации:

Метод последовательной релаксации для системы с пятидиа-гональной матрицей:

Время выполнения итерации в этом случае

.

Эти соотношения позволяют вычислить по столбцам элементы матрицы .

Эффективность такого метода достигается на этапе разложения матрицы, т. к. необходимо вычислить в этом случае только матрицу . Арифметические затраты в методе Холесского составляют длинных операций и n операций извлечения квадратного корня.

Существует другой вариант разложения симметрической положительно определенной матрицы, в котором удается избежать операций извлечения квадратного корня. В этом варианте вводится новая матрица по правилу

,

причем - матрица, вычисленная ранее по схеме Холесского, - диагональная матрица с элементами матрицы . Матрица существует и является нижней треугольной с единичной диагональю. В этом случае

,

Где

.

Расчетные соотношения для элементов матриц и можно получить, как и прежде, привлекая правило перемножения матриц

,

из которого следует, что

т. к. матрица имеет единичную диагональ.

Такой алгоритм потребует вдвое большего числа перемножений, чем схема Холесского. Однако, если ввести замену переменных

,

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

Здесь сначала вычисляют вспомогательные величины , а затем их используют для определения искомых величин и . Количество умножений при такой организации алгоритма составляет приблизительно и не содержит операции извлечения квадратного корня.

Лекция 5. Модификации метода сопряженных градиентов

ЦЕЛЬ ЛЕКЦИИ: Исследовать основные свойства метода сопряженных градиентов (СГ); построить на этой основе экономичную СГ-схему; показать, что переобусловливание матрицы решаемой системы позволяет повысить эффективность метода; получить вычислительную схему переобусловленного метода сопряженных градиентов.

5.1 Основные свойства метода сопряженных градиентов

Сопряженность. Свойство сопряженности означает, что

.

Это свойство в СГ-алгоритме обеспечивается выбором параметра .

Линейная независимость. Векторы , построенные по СГ-алгоритму, линейно независимы. Покажем это. Предположим обратное. Пусть один из векторов является линейной комбинацией остальных, т. е.

.

Умножим левую и правую части этого равенства слева сначала на A, затем на :

.

На основании свойства сопряженности . Отсюда следует ра-венство , что неверно, т. к. матрица A положительно опре-делена, т. е. для произвольного вектора выполняется неравен-ство . Значит, векторы линейно независимы.

Ортогональность невязки сопряженным градиентам:

.

В этом легко убедиться, если учесть, что

;

.

Умножим последнее равенство слева на и преобразуем его правую часть, привлекая свойство сопряженности ( для ):

.

Свойство доказано.

Ортогональность невязок:

.

Учтем, что

и, следовательно,

.

По этой причине

,

так как

и для .

5.2 Экономичная форма метода сопряженных градиентов

Преобразуем расчетные соотношения для и, используя перечисленные свойства ортогональности:

,

так как

и ;

вследствие того, что

В итоге метод сопряженных градиентов принимает экономичную форму:

1. Вычислить

.

2. Вычислить в цикле (k=0,1,2,…):

2.1. ; (- запомнить);

2.2. ;

2.3. ;

2.4. ;

2.5. .

Вычислительные затраты на итерации такой экономичной реализации метода сопряженных градиентов при решении СЛАУ с пятидиагональной матрицей сокращаются на время выполнения операций умножения и сложения:

.

Переобусловленный метод сопряженных градиентов.

Изменение ошибки на итерациях метода сопряженных градиентов характеризуется неравенством

,

где , - исходная погрешность и ее величина на i-й итерации; - число обусловленности матрицы решаемой СЛАУ. Чем ближе к единице, тем выше скорость сходимости метода. Ключ к построению более эффективного метода сопряженных градиентов сводится к следующему. Решаемая СЛАУ предварительно преобразуется таким образом, что матрица результирующей системы, сохраняя свойство симметричности и положительной определенности, характеризуется меньшим числом обусловленности. Такое преобразование называют переобусловливанием системы. К переобусловленной системе применяется стандартный метод сопряженных градиентов.

Каким образом можно осуществить переобусловливание системы? Введем симметрическую положительно определенную матрицу H такую, что, во-первых, , где - число обусловленности заключенной в скобках матрицы, во-вторых, решение линейной системы не требует существенных вычислительных затрат. Назовем матрицу переобусловливающей матрицей. Так как H - симметрическая матрица, то ее можно представить в виде , а переобусловливающую матрицу записать как

.

Формально переобусловливание системы с привлечением матрицы H осуществляется следующим образом:

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

причем

; ; ,

к которой и применяется метод сопряженных градиентов. Матрица симметрична, т. к. она получена из A с помощью преобразования конгруэнтности:, где - неособенная матрица.

Применим метод сопряженных градиентов к преобразованной системе:

1. Вычислить

.

2. Вычислить в цикле (k=0,1,2,…):

2.1.;

2.2.;

2.3.;

2.4.;

2.5..

Вернемся в этом алгоритме к переменной . Учтем, что. Тогда соотношение 2.2 алгоритма запишется следующим образом:

.

Если ввести вектор по правилу

,

то это соотношение переписывается в виде

,

т. е. становится итерационным правилом относительно x.

Преобразуем соотношение для :

.

Обозначим через невязку исходной системы. Тогда . Для произвольной k-й итерации можем записать:

.

В свою очередь

.

Соотношение для расчета невязки:

,

или

.

Выражение для :

.

Видоизменим правило построения вектора сопряженного градиента:

,

или

.

В итоге пришли к следующему алгоритму:

1. Построить матрицу .

2. Вычислить .

3. Вычислить в цикле (k=0,1,2,…):

3.1.;

3.2.;

3.3.;

3.4.;

3.5..

Реализация переобусловленного метода сопряженных градиентов предполагает построение матрицы . К настоящему времени разработан ряд способов построения матрицы H. Наиболее широко распространен на практике способ неполной факторизации, когда структура матрицы L совпадает со структурой нижней треугольной части матрицы A (либо близка к ней).

В частности, при решении систем линейных алгебраических уравнений с пятидиагональной матрицей в качестве матрицы H может использоваться матрица расщепления итерационного метода неполной факторизации.

5.3 Обобщенный метод сопряженных градиентов

Во многих практически важных приложениях возникает необходимость решения СЛАУ

,

матрица коэффициентов которой не является симметричной положительно определенной, т. е. и в зависимости от может быть как положительным, так и отрицательным числом. Применение метода сопряженных градиентов в таких случаях становится невозможным, так как функционал не является квадратичным. Тем не менее для решения подобных СЛАУ существует целый ряд методов, подобных методу сопряженных градиентов, которые в условиях точной арифметики гарантируют получение решения не более чем за n итерационных шагов.

Рассмотрим один из подходов к построению обобщенного метода сопряженных градиентов. Обратимся к минимизации следующей квадратичной функции

.

Из конечномерности пространства следует, что существует, доставляющее минимум данной функции, т. е.

.

Потребуем, чтобы матрица была положительно определена.

Будем строить итерационные приближения к следующим образом. Для заданного k определим k+1 векторов так, чтобы векторы были линейно независимы. Вдоль направлений определим k+1 параметр . Будем строить по итерационному правилу

,

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

.

Рассмотрим некоторое . Тогда условие записывается как . Преобразуем последнее равенство:

то есть

.

Перед тем как определить , отметим попутно, что

,

или окончательно

.

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

Соотношение для можно переписать как систему линей-ных уравнений

,

где

Будем строить векторы направлений в виде

,

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

.

Аналогично тому как это было сделано для традиционного метода сопряженных градиентов, умножим соотношение для вектора сопряженного направления слева сначала на A, а затем на:

.

Для левая часть здесь равна нулю. Справа под знаком суммы отличным от нуля будет лишь одно слагаемое, соответствующее . Тогда

.

Вернемся к уравнению для . В силу диагональности из всех отличным от нуля является коэффициент, соответствующий , а именно

.

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

Суммируя изложенное выше, запишем алгоритм обобщенного метода сопряженных градиентов:

Вычислить:

.

Вычислить в цикле ():

2.1.;

2.2. ;

2.3. ;

2.4. ;

2.5. .

В таком виде алгоритм соответствует обобщенному методу сопряженных невязок. К его достоинствам относятся, во-первых, свойство получения решения за конечное (не более n) число итерационных шагов, и, во-вторых, монотонное уменьшение нормы невязки с ростом . Данному методу присущи следующие недостатки. Во-первых, при расчете приходится использовать все k предыдущих векторов направлений, что при больших k делает метод трудно применимым в силу ограниченности объема оперативной памяти компьютера. Во-вторых, при расчете используются k предыдущих векторов.

На практике обобщенный метод сопряженных градиентов используется в одной из двух модификаций.

Задается целое число и выполняется s итераций. Затем вычисляется . Полагается и итерационный процесс перезапускается.

Одновременно хранится s последних векторов. При этом после расчета очередного последнее исключается. Нижний предел в сумме при расчете равен k-s+1.

...