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

2.4.2. Модификация метода Гаусса для случая линейных систем с трехдиагональными матрицами — метод прогонки

В приложениях часто возникают системы линейных уравнений с матрицами специального вида. В дальнейшем, например, при интерполяции функции сплайнами, при решении задачи Штурма-Лиувилля, при численном решении уравнений теплопроводности, уравнений Лапласа и Пуассона придется иметь дело с системами, матрицы которых имеют ненулевые элементы лишь на главной диагонали и еще на двух диагоналях — одной под главной, одной над главной. Такие матрицы называются трехдиагональными. Для решения систем с трехдиагональными матрицами существует экономичный (требующий малого количества памяти и арифметических действий) вариант метода Гаусса — прогонка. В англоязычной литературе метод прогонки называется алгоритмом Томаса. Подробнее о свойствах метода прогонки речь пойдет в лекции 10 в связи с решением разностными методами задачи Штурма-Лиувилля.

2.4.3. LU-разложение

Среди прямых методов численного решения СЛАУ широко используется также LU-разложение матрицы $$\mathbf{A}$$ и метод Холецкого (или метод квадратного корня).

Если матрица $$\mathbf{A}$$ представима в виде произведений матриц $$\mathbf{LU}$$, то СЛАУ может быть представлена в виде

Перепишем (2.13), вводя вспомогательный вектор $$\mathbf{v}$$, в следующем виде

Решение СЛАУ свелось к последовательному решению двух систем с треугольными матрицами. Первый этап решения системы $$\mathbf{Lv} = \mathbf{f}$$:

откуда можно вычислить все vk последовательно по формулам

Далее рассмотрим систему $$\mathbf{Uu} = \mathbf{v}$$ или

решение которой находится в обратном порядке, т.е. при k = n - 1, ..., 1 по очевидным формулам $$u_k = d_{kk}^{- 1}(v_k - \sum\limits_{j = k + 1}^n{d_{kj}u_j )}.$$ Условия существования такого разложения даются следующей теоремой [2.5] (без доказательства).

Теорема. Если все главные миноры квадратной матрицы $$\mathbf{A}$$ отличны от нуля, то существуют единственные нижняя и верхняя треугольные матрицы $$\mathbf{L} = l_{ij}$$ и $$\mathbf{U} = d_{ij}$$ такие, что $$\mathbf{A}= \mathbf{LU}.$$ При этом все диагональные коэффициенты матрицы $$\mathbf{L}$$ фиксированы и равны единице.

Опишем алгоритм нахождения элементов lijdij матриц $$\mathbf{L}, \mathbf{U}.$$ Выписав равенство $$\mathbf{A} = \mathbf{LU}$$ в компонентах, получим

Выполнив умножение матриц, приходим к системе линейных уравнений размером n x n:

относительно неизвестных d11, d12, ..., d1n, l21, d22, ..., d2n, ln1, ln2, ..., dnn.

Специфика этой системы позволяет решить ее последовательно. Из первой строки находим d1j= a1j(j = 1, ..., n).

Из уравнений, входящих в первый столбец приведенной выше системы, находим li1= ai1/d11, i = 1, ...dots, n. Теперь можно из уравнений второй строки найти d2j= a2j- l21d1j, j = 2, ..., n, а из уравнений, входящих во второй столбец, получим $$l_{i2}= d_{22}^{-1}(a_{i2}- l_{i1}d_{12}), i = 2, \ldots, n$$ и так далее. Последним вычисляется элемент

Можно выписать общий вид этих формул:

Приведение матриц к треугольному виду аналогично приведению матрицы в методе Гаусса и также требует количества арифметических действий порядка O(n3), точнее, $$\approx 2n^{3}.$$

2.4.4. Метод Холецкого (метод квадратного корня)

Пусть матрица рассматриваемой линейной системы $$\mathbf{A}$$ — симметричная, т.е. aij = aji, положительная матрица. Тогда она представима в виде $$\mathbf{A} = \mathbf{LL}^T$$, где

Далее, как и в случае LU-разложения, решение СЛАУ $$\mathbf{Au} = \mathbf{f}$$ сводится к последовательному решению двух линейных систем с треугольными матрицами $$\mathbf{Lv} = \mathbf{f}$$, $${\mathbf{L}}^T{\mathbf{u}} = \mathbf{v}$$, для решения которых требуется примерно 2n2 арифметических действий.

Первая из этих линейных систем

она легко решается. Для решения получаем очевидные формулы

Вторая система уравнений есть

l11u1  + l12u2  +  ...  + l1nun  = v1 , l22u2  +  ...  + l2nun  = v2, ... lnnun  = vn.

Из нее находим значения переменных ui в обратном порядке по формуле

Определенной опасностью при реализации этого метода являются возможная близость к нулю lii и отрицательность подкоренных выражений при вычислении lii (последнего не должно быть при симметричной положительной матрице $$\mathbf{A}$$ )

Элементы матрицы $$\mathbf{L}$$ находим из уравнения $${\mathbf{LL}}^T = \mathbf{A}$$, приравнивая соответствующие элементы матриц $${\mathbf{LL}}^T$$ и $$\mathbf{A}.$$ В результате получим систему уравнений

Решение этой системы легко находится:

Метод также называется методом квадратного корня.

Внимание! Не следует путать матрицу (оператор) $$\mathbf{L}$$ с оператором $$\mathbf{A}^{1/2}$$ — квадратным корнем из самосопряженного положительного оператора.

2.5. Итерационные методы решения СЛАУ

2.5.1. Метод простой итерации

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

Проведем несколько равносильных преобразований. Умножим обе части системы на один и тот же скалярный множитель $$\tau,$$ затем прибавим к правой и левой частям системы вектор $$\mathbf{u}.$$ Систему уравнений можно теперь записать в виде, удобном для итераций:

где $$\mathbf{B} = \mathbf{E} - \tau \mathbf{A},\quad \mathbf{F} = \tau \mathbf{f}.$$

Теперь построим последовательность приближений к решению системы. Выберем произвольный вектор $${\mathbf{u}}_0$$ — начальное приближение к решению. Чаще всего его просто полагают нулевым вектором. Скорее всего, начальное приближение не удовлетворяет (2.15) и, следовательно, исходной системе. При подстановке его в исходное уравнение возникает невязка $${\mathbf{r}}_0 = \mathbf{f} - \mathbf{A}{\mathbf{u}}_0.$$ Вычислив невязку, с помощью (2.15) можно уточнить приближение к решению, считая, что

По первому приближению снова вычисляется невязка, процесс продолжается. В ходе итерации получаем $${\mathbf{u}}_{k + 1} = {\mathbf{u}}_k + \tau{\mathbf{r}}_k, {\mathbf{r}}_k = \mathbf{f} - {\mathbf{Au}}_k.$$ Эквивалентная формулировка метода, называемого методом простых итераций, заключается в следующем. Решение (2.15) находится как предел последовательности $$\{{\mathbf{u}}_0, {\mathbf{u}}_1, {\mathbf{u}}_2, \ldots\}$$ приближений, члены которой связаны рекуррентным соотношением (оно эквивалентно приведенному выше, из записи исключен вектор невязки):

$${\mathbf{u}}_0 = 0$$ (или любому произвольному вектору). Если предел такой последовательности существует, то говорят о сходимости итерационного процесса к решению СЛАУ.

Существуют другие формы записи метода итераций, например

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

При $${\mathbf{D}}_k = {\mathbf{E}}$$, $$\tau _{k} = \tau$$ последняя формула соответствует однопараметрическому итерационному процессу — рассмотренному выше методу простых итераций. При $${\mathbf{D}}_k = \mathbf{E}$$, $$\tau_k= \left\{{\tau_k, k= 1, \ldots , n}\right\}$$ — n -шаговому явному итерационному процессу, при $${\mathbf{D}}_k = \mathbf{D^{\prime}}$$, $$\tau _{k} = 1$$ — методу простой итерации без итерационного параметра. В случае, когда $$\mathbf{D} \ne \mathbf{E}$$, итерационный метод называется неявным — для вычисления следующего приближения к решению придется решать (как правило, более простую, чем исходную) систему линейных уравнений.

Теорема (достаточное условие сходимости метода простой итерации ). Итерационный процесс (2.16) сходится к решению $$\mathbf{U}$$ СЛАУ $$\mathbf{Au} = \mathbf{F}$$ со скоростью геометрической прогрессии при выполнении условия: $$\|{\mathbf{B}}\| \le q < 1.$$

Доказательство.

Пусть $$\mathbf{U}$$ — точное решение системы (2). Вычитая из (2.16)-(2.15), получим $${\mathbf{u}}_k - {\mathbf{U}} = \mathbf{B}({\mathbf{u}}_{k - 1} - \mathbf{U})$$, или, обозначив погрешность $${\mathbf{\varepsilon }}_k = {\mathbf{u}}_k - {\mathbf{U}}$$, получим для эволюции погрешности уравнение $${\mathbf{\varepsilon}}_k = {\mathbf{B\varepsilon}}_{k - 1}.$$ Справедлива цепочка неравенств: $$\|{{\mathbf{u}}_k - {\mathbf{U}}}\| = \|{{\mathbf{\varepsilon}}_k}\| \le \|{\mathbf{B}}\| \cdot \|{{\mathbf{\varepsilon}}_{k - 1}}\| \le q\|{{\mathbf{\varepsilon}}_{k - 1}}\| \le \ldots \le q^k\|{{\mathbf{\varepsilon}}_0}\| = q^k\|{{\mathbf{u}}_0 - {\mathbf{U}}}\|$$, где $$0 < q \le \|{\mathbf{B}}\|.$$

Отсюда следует, что при $$q < 1 \lim\limits_{k \to \infty}{\mathbf{u}}_k = \mathbf{U}.$$

Из неравенства $$\|{{\mathbf{\varepsilon}}_k}\| \le q^k\|{{\mathbf{\varepsilon}}_0}\|$$ можно получить оценку количества итераций,необходимых для достижения точности $$\varepsilon,$$ т.е. для выполнения условия $$\|{{\mathbf{u}}_k - \mathbf{U}}\| = \|{{\mathbf{\varepsilon}}_k}\| \le \varepsilon.$$ Эта оценка имеет вид $$$ k \ge \left({\ln{\frac{\mathbf{\varepsilon}}{\|{{\mathbf{\varepsilon}}_0}\|}}}\right)/\ln{q} $.$$

Теорема (критерий сходимости метода простой итерации (без доказательства)). Пусть СЛАУ (2.15) имеет единственное решение. Тогда для сходимости итерационного процесса (2.16) необходимо и достаточно, чтобы все собственные значения матрицы $${\mathbf{B}$$ по абсолютной величине были меньше единицы.

Сравним по количеству арифметических действий прямые и итерационные методы. Метод Гаусса без выбора главного элемента при $$n \gg 1$$ требует $$$ \approx (\frac{2}{3}n^3) $$$ арифметических действий; метод простой итерации (2.16) $$\approx (2n^2 \cdot I)$$, где i — число приближений, необходимое для достижения заданной точности. Значит, при I < n/3 метод итераций становится предпочтительнее. В реальных задачах, в основном, $$I \ll n.$$ Кроме того, итерационные методы можно делать более эффективными, изменяя итерационные параметры. В ряде случаев итерационные методы оказываются более устойчивыми по отношению к накоплению ошибок округления, чем прямые.

2.5.2. Влияние ошибок округления на результат численного решения

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

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

где $${\mathbf{\delta }}_k$$ — суммарная погрешность округления. Норму разности между реальным и идеальным (т.е. в отсутствии ошибки округления) результатами расчетов получим, вычитая (2.19) из (2.20). Учтем, что $$\|\mathbf{B}\| < q < 1$$,

Так как начальное приближение задано точно $$\|{{\mathbf{u}}_0^M - {\mathbf{u}}_0}\| = 0.$$ Обозначим $$\delta = \max\limits_i \|{{\mathbf{\delta }}_i}\|$$ и вычислим сумму членов геометрической прогрессии. Получим $$$ \|{{\mathbf{u}}_k^M - {\mathbf{u}}_k}\| \le \delta \frac{q^k - 1}{q - 1} \le \frac{\delta }{1 - q} $$$, то есть погрешность, вносимая в решение из-за конечной разрядности мантиссы, не зависит от количества итераций. Этот результат является характеристикой устойчивости рассматриваемого вычислительного процесса.

2.5.3. Методы Якоби, Зейделя, верхней релаксации

Представим матрицу $$\mathbf{A}$$ в виде

где $$\mathbf{L}$$ и $$\mathbf{U}$$ — нижняя и верхняя треугольные матрицы с нулевыми элементами на главной диагонали, $$\mathbf{D}$$ — диагональная матрица. Рассматриваемая СЛАУ может быть переписана в следующем эквивалентном виде:

Построим два итерационных метода $$\mathbf{Lu}_k + \mathbf{Du}_{k+1} + \mathbf{Uu}_k = \mathbf{f} $$

и

или, соответственно,

и

Очевидно, что эти формулы описывают итерационные процессы вида (2.16), если положить в (2.22)

или

Эти итерационные процессы называются методами Якоби и Зейделя . Представим их в компонентной записи. Метод Якоби будет иметь вид (перенесем итерационный индекс k вверх):

Метод Зейделя можно представить следующим образом:

Эти формулы легко выводятся, если учесть, что элементами матрицы D -1 являются $$d_{ii} = a_{ii}^{- 1}.$$

Теорема (достаточное условие сходимости метода Якоби ) Итерационный метод Якоби сходится к решению соответствующей СЛАУ, если выполнено условие диагонального преобладания

Доказательство.

Выполненные условия (2.24) означает, что в любой строке матрицы перехода

сумма модулей элементов меньше единицы. В этом случае по крайней мере одна из норм матрицы $$\mathbf{B}$$ меньше единицы. Тогда выполняется достаточное условие сходимости метода простых итераций.

Теорема (критерий сходимости итерационного метода Якоби ). Для сходимости итерационного метода Якоби необходимо и достаточно, чтобы все корни уравнения

по модулю не превосходили единицы.

Доказательство.

Легко проверить, что в силу диагональности $$\mathbf{D}$$ имеет место

Собственными значениями матрицы $$\mathbf{B} = -\mathbf{D}^{- 1}(\mathbf{L} + \mathbf{U})$$ являются корни уравнения

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

Аналогичную теорему можно доказать и для метода Зейделя, однако матрица в этой теореме будет иметь другой вид:

Теорема (достаточное условие сходимости метода Зейделя (без доказательства)). Пусть $$\mathbf{A}$$ — вещественная, симметричная, положительно определенная матрица. В этом случае итерационный метод Зейделя сходится.

Доказательство этой теоремы сводится к проверке того, что выполнение условий теоремы для матрицы $$\mathbf{A}= \mathbf{L} + \mathbf{D} + {\mathbf{L}}^{T}$$ влечет выполнение условия сходимости итерационного метода с матрицей перехода — $${(\mathbf{L}+\mathbf{D})}^{-1}\mathbf{L}^T.$$ СЛАУ с вещественной матрицей $$\mathbf{A}$$ такой, что $$\det \mathbf{A} \ne 0$$ может быть симметризована умножением на матрицу $${\mathbf{A}}^T$$:

(симметризация Гаусса).

Развитием метода Зейделя является метод релаксации. В этом методе вводится итерационный параметр $$\tau,$$ называемый параметром релаксации. Представим метод релаксации в матричной форме:

Выбирая $$\tau,$$ можно существенно изменять скорость сходимости итерационного метода. Выразим $${\mathbf{u}}_{k+1}$$

В общем случае задача вычисления $$\tau _{опт}$$ (оптимального итерационного параметра) не решена, однако известно, что $$1 < \tau _{опт} < 2.$$ В этом случае итерационный метод называется методом последовательной верхней релаксации или SОR — Successive Over Relaxation. Иногда встречается термин "сверхрелаксация" при $$1 < \tau _{опт} < 2.$$ При $$0 < \tau < 1$$ имеем метод нижней релаксации.

2.6. Вариационные итерационные методы

2.6.1. Связь между вариационной задачей и задачей решения СЛАУ

Пусть $$\mathbf{u} \in L^n $$, где Ln есть n -мерное евклидово пространство. Рассмотрим квадратичный функционал от $$\mathbf{u}$$, называемый функционалом энергии:

где $$\mathbf{A}$$ — линейный оператор, $$\mathbf{f} \in L^n$$, c — константа. Этот функционал совпадает с квадратичным функционалом $$\Phi (\mathbf{u}) = (\mathbf{A}^*\mathbf{u,u}) - 2(\mathbf{f,u}) + c$$, где $$\mathbf{A}^*$$ — сопряженный к $$\mathbf{A}$$ оператор. Действительно, $$(\mathbf{Au,u}) \equiv (\mathbf{u},\mathbf{A}^*\mathbf{u})$$ по определению сопряженного оператора и $$(\mathbf{u},\mathbf{A}^*\mathbf{u}) = (\mathbf{A}^*\mathbf{u},\mathbf{u})$$ в силу коммутативности скалярного произведения. Тогда

Без ограничения общности предположим, что оператор $$\mathbf{A}$$ — самосопряженный, $$\mathbf{A} = \mathbf{A}^*.$$ В противном случае будем рассматривать задачу с оператором $$\frac{1}{2}(\mathbf{A}+{\mathbf{A}^*}) $$ при решении вариационной задачи.

Будем также считать, что $$\mathbf{A}$$ — положительный оператор, т.е. $$\mathbf{A} > 0$$, это означает, что для любого ненулевого вектора $$\mathbf{u}$$ выполнено $${(\mathbf{Au}, \mathbf{u}) > 0}.$$

Поставим задачу об отыскании элемента $$\mathbf{v}$$, придающего наименьшее значение функционалу $$\Phi (\mathbf{u})$$:

Теорема. Пусть $$\mathbf{A} = {\mathbf{A}*} > 0.$$ В этом случае существует единственный элемент $${\mathbf{v}} \in L^n $$, придающий наименьшее значение квадратичному функционалу $$\Phi (\mathbf{u}) = (\mathbf{Au,u}) - (2\mathbf{f,u}) + c$$, являющийся решением СЛАУ $$\mathbf{Au}= \mathbf{f}.$$

Доказательство.

СЛАУ $$\mathbf{Au}= \mathbf{f}$$ имеет единственное решение $$\mathbf{v}$$, поскольку $$\mathbf{A}$$ является невырожденным оператором в силу его положительной определенности. Покажем, что в этом случае при $${\mathbf{Av}} - \mathbf{f} = 0$$ для любого вектора $$\Delta$$ имеет место $$\Phi ({\mathbf{v}} + {\mathbf{\Delta }}) > \Phi ({\mathbf{v}})$$, т.е. при $$\mathbf{u} = \mathbf{v}$$ достигается минимум квадратичного функционала $$\Phi (\mathbf{u}).$$

Действительно,

т.е. при $$\mathbf{Av} = \mathbf{f}$$ и любом $$\Delta$$ имеет место $$\min\limits_\mathbf{u}\Phi (\mathbf{u}).$$ Докажем, что верно и обратное утверждение. Если элемент доставляет минимальное значение функционалу энергии, то он является решением системы линейных уравнений $$\mathbf{Av} = \mathbf{f}.$$ Из курса математического анализа известно, что в точке минимума должно выполняться условие $$grad \Phi (\mathbf{u}) = 0,\quad \mathbf{A} > 0.$$ Вычисляя градиент, приходим к условию минимума функционала $$grad \Phi (\mathbf{u}) = 2\mathbf{Au}- 2\mathbf{f}= 0.$$ Таким образом установлена эквивалентность вариационной задачи (отыскание элемента, придающего минимум $$\Phi (\mathbf{u})$$ ) и задачи о нахождении решения СЛАУ.

Заметим, что СЛАУ с самосопряженным и положительно определенным оператором $$\mathbf{A}$$ представляют собой важный класс задач в математической физике, в частности, они возникают при решении краевых задач для эллиптических уравнений. При необходимости можно произвести симметризацию по Гауссу исходной системы.

2.6.2. Методы градиентного и наискорейшего спуска

Метод градиентного спуска состоит в нахождении следующего приближения в итерационном процессе из предыдущего, путем смещения в направлении градиента функционала

где $$\mathbf{A}$$ — положительно определенная симметричная матрица; $${\alpha}_k$$ — параметр, определяемый из заданных условий; например, из условия минимума величины

В этом случае итерационный метод называется методом наискорейшего спуска. Так как $$grad \Phi (\mathbf{u}) = 2(\mathbf{Au}- \mathbf{f})$$, то (2.26) приобретает вид

что соответствует записи итерационного метода в форме {2.17). Здесь $$\tau _{k}$$ является итерационным параметром, который в методе наискорейшего спуска определяется из условия минимума функции $$\Phi (\tau_k, {\mathbf{u}}_{k + 1})$$ по $$\tau _{k}.$$ Найдем условие этого минимума:

Здесь учтено соотношение: $$({\mathbf{Au^{\prime}}},\mathbf{u}) = ({\mathbf{u^{\prime}}},\mathbf{A}^*\mathbf{u}) = (\mathbf{Au}{\mathbf{,u^{\prime}}})$$, поскольку $$\mathbf{A} = \mathbf{A}^*$$ и $$(\mathbf{v,Aw}) = (\mathbf{Aw,v})$$ в силу самосопряженности оператора $$\mathbf{A}.$$ Подставим в последние равенства $${\mathbf{u}}_{k+1}$$ из (2.27), получим $$(\mathbf{Au}_k - \mathbf{f}- \tau_k\mathbf{A}(\mathbf{Au}_k - \mathbf{f}),\mathbf{Au}_k - \mathbf{f}) = 0$$, откуда следует

$$\tau_k = (\mathbf{Au}_k - \mathbf{f},\mathbf{Au}_k - \mathbf{f})/(\mathbf{A}(\mathbf{Au}_k - \mathbf{f}),\mathbf{Au}_k - \mathbf{f})$$, или $$$ \tau_k = \frac{({\mathbf{r}}_k,{\mathbf{r}}_k)}{({\mathbf{Ar}}_k,{\mathbf{r}}_k)}$$, где $${\mathbf{r}}_k = \mathbf{Au}_k - \mathbf{f} $.$$

Вектор $${\mathbf{r}}_k$$ называют вектором невязки.

2.6.3. Метод минимальных невязок

Этот итерационный метод определяется следующим образом. Пусть $${\mathbf{u}}_{k+ 1} = {\mathbf{u}}_k - \tau_k{\mathbf{r}}_k$$, как и ранее, $${\mathbf{r}}_k = \mathbf{Au}_k - \mathbf{f}.$$ Итерационный параметр $$\tau _{k}$$ на каждой итерации выбирается так, чтобы минимизировать евклидову норму невязки $${\mathbf{r}_{k+1}}.$$ Заметим, что итерационный процесс $${\mathbf{u}}_{n + 1} = {\mathbf{u}}_n + \tau_n{\mathbf{r}}_{n}$$ может быть представлен в равносильном виде в терминах невязки $${\mathbf{r}}_{n + 1} = {\mathbf{r}}_{n} + \tau_n{\mathbf{Ar}}_n.$$ Тогда для квадрата евклидовой (третьей) нормы невязки получаем условие

Для отыскания минимума невязки на следующей итерации приравняем нулю производную последнего выражения по итерационному параметру $$\tau _{k}.$$ Получим равенство

Из последнего соотношения находим значение итерационного параметра $$$ \tau_k= \frac{({\mathbf{Ar}}_k,{\mathbf{r}}_k)}{({\mathbf{Ar}}_k,{\mathbf{Ar}}_k)}. $$$

2.6.4. Метод сопряженных градиентов

Этот метод применяется для решения систем уравнений с самосопряженной положительной матрицей $$\mathbf{A} = \mathbf{A}^* > 0.$$ Оптимизируем градиентный метод, выбирая параметры $$\tau$$ таким образом, чтобы на последующем шаге невязка была ортогональна всем предыдущим. На первом шаге невязку ищем аналогично методу наискорейшего спуска. Получим невязки, образующие ортогональный базис. На последнем шаге невязка становится равна нулю, так как пространство конечномерно, и единственный элемент, ортогональный всем базисным векторам конечномерного пространства — нулевой. Получаем точное решение за конечное число шагов ( прямой метод ). Однако этот метод работает не всегда, так как при плохой обусловленности матрицы он становится вычислительно неустойчивым.

Идея метода состоит в следующем. Выбираем произвольное начальное приближение и вычисляем по нему вектор невязки $${\mathbf{r}}^0 = \mathbf{Au}^0 - \mathbf{f}$$, тогда первое приближение $$\mathbf{u}^1 = \mathbf{u}^0 + \tau_0{\mathbf{r}}^0 .$$

Из условия ортогональности невязок на двух первых шагах находим значение итерационного параметра $$$ \tau_0 = - \frac{({\mathbf{r}}^0,{\mathbf{r}}^0 )}{({\mathbf{Ar}}^0 ,{\mathbf{r}}^0 )} $.$$

Построим такое приближение, чтобы учитывались две предыдущие — трехслойный итерационный метод. Фактически, при построении его применяется процесс ортогонализации Грамма-Шмидта. Если $$\left({{\mathbf{Ar}}^{n - 2},{\mathbf{r}}^{n - 1}}\right) = 0$$, то $${\mathbf{r}}^{n - 2},{\mathbf{r}}^{n - 1}$$ — $$\mathbf{A}$$ - сопряженные невязки.

Имеем метод сопряженных градиентов.

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

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

где

В вычислительной практике этот метод используется при умеренном числе обусловленности, больших n и неизвестных границах спектра матрицы $$\mathbf{A}$$, как итерационный метод, поскольку k0 обычно достаточно большое число. Подробнее о методах сопряженных градиентов можно прочитать в [2.2].

2.7. О спектральных задачах

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

Рассмотрим два алгоритма для самосопряженных матриц. Первый — степенной алгоритм, для вычисления наибольшего по абсолютной величине собственного числа. Выбираем произвольный ненулевой вектор $${\mathbf{u}}_0$$ и строим последовательность векторов

Легко показать, что выражение

приближает максимальное по абсолютной величине собственное значение с точностью $$O(\lambda_{N}/\lambda_{N- 1})^k.$$ Здесь $$\lambda _{N}/\lambda _{N - 1}$$ — отношение самого большого по модулю собственного числа матрицы к следующему по абсолютной величине.

Для решения полной самосопряженной проблемы собственных значений применяется метод вращений.

Определение собственных значений самосопряженной матрицы $$\mathbf{A}$$ эквивалентно отысканию такой ортогональной матрицы $$\mathbf{T}$$, что

матрица $$\mathbf{\Lambda }$$ — диагональная. Среди всех ортогональных преобразований данное минимизирует сумму квадратов внедиагональных элементов исходной матрицы. Построим итерационный метод, минимизирующий эту сумму на каждой итерации. Пусть каждое преобразование подобия на каждой итерации содержит лишь одну матрицу вращения $${\mathbf{\hat A}} = {\mathbf{T^{\prime}}}_{ij}{\mathbf{AT}}_{ij}$$, где матрица $$\mathbf{T}_{ij}$$ есть матрица поворота в плоскости $${\mathbf{u}}_i {\mathbf{u}}_j$$ на угол $$\alpha.$$ Эта матрица отличается от матрицы $$\mathbf{A}$$ только двумя строками и двумя столбцами (с номерами i и j ). Так как евклидова норма матрицы не изменяется при ортогональных преобразованиях, то легко получить соотношение между суммами квадратов внедиагональных элементов старой и новой матриц:

Очевидны условия минимизации суммы в левой части последнего равенства. Следует на текущей итерации выбирать индексы так, чтобы выполнялось условие $$|{a_{ij}}| = \max\limits_{k \ne l}|{a_{kl}}|$$, а угол поворота выбирается из условия $$0 = {\left({(a_{jj}} - a_{ii})\sin 2\alpha + 2a_{ij}\cos 2\alpha}\right)}^2 .$$ Тогда он удовлетворяет условию $$$ \tg 2\alpha = \frac{2a_{ij}}{a_{ii} - a_{jj}},\quad |{\alpha }| \le \frac{\pi}{4} $.$$

Независимо от наличия кратных собственных значений метод вращений обладает квадратичной сходимостью.

Выбор максимального по модулю внедиагонального элемента — затратная операция, поэтому часто реализуется метод вращений с барьерами. Его идея состоит в следующем. При переборе внедиагональных значений вращение производится тогда, когда значение элемента по абсолютной величине превосходит некоторую величину (барьер). Если все элементы меньше барьера, его значение уменьшается, например, на порядок, и снова начинается циклический перебор внедиагональных элементов. Подробнее о методе вращений смотри в [2.1].

Другие алгоритмы решения спектральных задач описаны в специальной литературе [2.1, 2.3, 2.10].

Задачи

Методом Гаусса решить систему линейных уравнений $$\mathbf{Ax} = \mathbf{f}$$, где

$$\mathbf{A} ={\left( \begin{array}{ccc} 5 0 1 \\ 2 6 {-2}\\ -3 2 10 \\ \end{array}\right)}, \mathbf{f} = {(11, 8, 6)}^{\mathbf{T}}$$

или

$$\begin{gather*} 5{x}_1 + 0 \cdot{x}_2 +{x}_3 = 11, \\ 2{x}_1 + 6 \cdot{x}_2 + 2{x}_3 = 8, \\ - 3{x}_1 + 2 \cdot{x}_2 + 10{x}_3 = 6. \end{gather*} $$

Решение.Расширенная матрица $$$ \tilde{\mathbf{A}}$$$ имеет вид

$$$\tilde{\mathbf{A}} = \left( \begin{array}{cccc} {{a}_{11}} {{a}_{12}} {{a}_{13}} {f}_1 \\ {{a}_{21}} {{a}_{22}} {{a}_{23}} {f}_2 \\ {{a}_{31}} {{a}_{32}} {{a}_{33}} {f_3}\\ \end{array}\right) = \left( \begin{array}{cccc} 5 0 1 11\\ 2 6 {-2} 8 \\ -3 2 10 6 \\ \end{array}\right) $$$

Разделив элементы первой строки на ведущий элемент a11= 5 получаем первую опорную строку (1, 0, 0.2, 2.2).

Далее умножим ее на a21= 2 и вычтем из второй строки, после чего умножим опорную строку на a31= 3 и вычтем из третьей. Получаем матрицу $$$ \tilde{\mathbf{A}}_1 $$$

$$$ \tilde{\mathbf{A}} = \left( \begin{array}{cccc} {{a}_{11}} {{a}_{12}} {{a}_{13}} {{f}_1}\\ 0 {\tilde{a}_{22}} {\tilde{a}_{23}} {\tilde{f}_2}\\ 0 {\tilde{a}_{32}} {\tilde{a}_{33}} {\tilde{f}_3}\\ \end{array}\right) = \left( \begin{array}{cccc} 1 0 0.2 2.2\\ 0 6 {-2.4} 3.6 \\ 0 2 10.6 12.6 \\ \end{array}\right) $$$

Вторая опорная строка — результат деления второй строки матрицы $$$ \tilde{\mathbf{A}}_1 $$$ на

$$a_{22}^1 = 6$$:

(0, 1, - 0.4, 0.6)

Матрицы $$$ \tilde{\mathbf{A}}_2 $$$ после умножения опорной строки на 2 и вычитания ее из третьей

$$$ \tilde{\mathbf{A}}= \left( \begin{array}{cccc} {{a}_{11}} {{a}_{12}} {{a}_{13}} {f}_1 \\ 0 {\tilde{a}_{22}} {\tilde{a}_{23}} {\tilde{f}_2}\\ 0 0 {\tilde{\tilde{a}}_{33}} {\tilde{\tilde{f}}_3}\\ \end{array}\right) = \left( \begin{array}{cccc} 1 0 0.2 2.2\\ 0 6 {-0.4} 0.6 \\ 0 0 11.4 11.4 \\ \end{array}\right) $$$

Третья опорная строка (результат деления третьей строки на 11,4 ) есть (0,0,1,1) .

Матрица $$$ \tilde{\mathbf{A}}_3 $$$ будет

$$$ \tilde{\mathbf{A}}_1 = \left( \begin{array}{cccc} 1 0 {0,2} {2,2}\\ 0 1 {- 0,4} {0,6}\\ 0 0 1 1 \\ \end{array}}\right). $$$

Обратный ход метода Гаусса.

Система уравнений с матрицей $$$ \tilde{\mathbf{A}}_3 $$$

$$\begin{gather*} {x}_1 + 0 \cdot{x}_2 + 0,2{x}_3 = 2,2 , \\ {x}_2 - 0,4{x}_3 = 0,6 , \\ {x}_3 = 1. \end{gather*} $$

Разрешая эту систему, начиная с последнего уравнения, получим

x3  = 1,  x2  = 1,  x1  = 2.

Определитель матрицы $$\det \mathbf{A}$$ можно вычислить как произведение ведущих элементов

$$\det \mathbf{A} = 5 \cdot 6 \cdot 11,4 = 342.$$

Показать, что решение системы линейных уравнений

$$\begin{gather*} {x}_1 + 2{x}_2 + 3{x}_3 + 4{x}_4 = 2, \\ {x}_1 + 3{x}_2 + {x}_3 + 2{x}_4 = - 1, \\ 2{x}_1 + 3{x}_2 + 8{x}_3 + 7{x}_4 = 10, \\ 2{x}_1 + 5{x}_2 + 3{x}_3 + 7{x}_4 = 3 \end{gather*}$$

методом Гаусса невозможно.

Решение. Преобразование матрицы рассматриваемой системы будет

$$\left( \begin{array}{ccccc} 1 2 3 4 2 \\ 1 3 1 2 {- 1}\\ 2 3 8 7 {10}\\ 2 5 3 7 3 \\ \end{array}\right) \to \left( \begin{array}{ccccc} 1 2 3 4 2 \\ 0 1 {- 2} {- 2} {- 3}\\ 0 {- 1} 2 {- 1} 6 \\ 0 1 {- 3} {- 1} 1 \\ \end{array}\right) \to \left( \begin{array}{ccccc} 1 2 3 4 2 \\ 0 1 {- 2} {- 2} {- 3}\\ 0 0 0 {- 3} {- 3}\\ 0 0 {- 1} 1 2 \\ \end{array}\right)$$

Поскольку $${a^2_{33}} = 0$$ то вычисление третьей опорной строки невозможно.

Показать, что

$$$ \mu (\mathbf{A}) \ge \frac{|{\max\limits_i \lambda_i(\mathbf{A})}|}{|{\min\limits_i \lambda_i (\mathbf{A})}|}. $$$

Рассмотреть случай симметричной матрицы $$\mathbf{A}.$$

Решение. Для собственного вектора $$\omega,$$ соответствующего наибольшему по модулю собственному значению матрицы, выполняется равенство $$\mathbf{A}\omega = \lambda \omega $$, откуда

$$\|{\mathbf{A\omega}}\| = \left|{\max\limits_i\lambda_i}\right|{\|\omega\|}.$$

Учитывая, что $$\|{\mathbf{A}\omega }\| \le \|\mathbf{A}\|{\|\omega\|}$$, получим $$\|\mathbf{A}\| \ge \left|{\max\limits_i \lambda_i (\mathbf{A})}\right|.$$ Для обратной матрицы $$\mathbf{A}^{-1}$$ максимальным по модулю является собственное число $$\min\limits_i \lambda_i^{-1}$$, откуда $$\left\|{\mathbf{A}^{- 1}}\right\| \ge {\left|{\min\limits_i\lambda_i(\mathbf{A})}\right|}^{- 1}.$$ Объединяя два последних неравенства, получим

$$$ \mu (\mathbf{A}) = {\left\|{\mathbf{A}^{- 1}}\right\|}\|\mathbf{A}\| \ge \frac{|{\max\limits_i \lambda_i (\mathbf{A})}|}{|{\min\limits_i \lambda_i (\mathbf{A})}|}. $$$

В случае симметричной матрицы $$\mathbf{A}(\mathbf{A}^* = \mathbf{A})$$ имеем

$$\|\mathbf{A}\|}_3 = \sqrt{\lambda_{\max}(\mathbf{A}^*\mathbf{A})} = \sqrt{\lambda_{\max}(\mathbf{A}^2 )}= \sqrt{\lambda_{\max}^2(\mathbf{A})} = |{\lambda_{\max}(\mathbf{A})}|$$,

т.к. из $${\mathbf{A\omega}}_i = \lambda_i\omega_i$$ следует $$({\mathbf{A\omega}}_i )^2 = \lambda_i^2\omega_i^2 .$$ Аналогично

$${\left\|{\mathbf{A}^{- 1}}\right\|}_3 = \sqrt{\lambda_{\max}[(\mathbf{A}^{- 1})^* (\mathbf{A}^{- 1})]} = \sqrt{\lambda_{\max}[(\mathbf{A}^{- 1})]^2} = \\ = \sqrt{\lambda_{\min}^{- 1}(\mathbf{A}^2 )} = |{\lambda_{\min}^{- 1}(\mathbf{A})}|.$$

Тогда $$\mu (\mathbf{A}) = {\|\mathbf{A}\|}_3{\left\|{\mathbf{A}^{- 1}}\right\|}_3 = |{\lambda_{\max}(\mathbf{A})}|/|{\lambda_{\min}(\mathbf{A})}|.$$

Найти число обусловленности матрицы $$\mathbf{A}$$, выразив его через число обусловленности матрицы $$\mathbf{B}$$, если $$\mathbf{A} = {\mathbf{B}}^*{\mathbf{B}} > 0.$$

Решение. Для самосопряженной положительной матрицы $$\mathbf{A}$$ имеем

$$$ \|\mathbf{A}\| = \sup\limits_{\|\mathbf{u}\| \ne 0}\frac{(\mathbf{u}, \mathbf{Au})}{(\mathbf{u}, \mathbf{u})}. $$$

Тогда

$$$ \mu (\mathbf{A}) = {\|\mathbf{A}\|}_3 \left\|{{\mathbf{A}}^{- 1}}\right\|}_3 = \sup\limits_{\|\mathbf{u}\| \ne 0}\frac{(\mathbf{u, Au})}{(\mathbf{u, u})} \sup\limits_{\|{u}\| \ne 0}\frac{(\mathbf{u}, \mathbf{A}^{- 1}\mathbf{u})}{(\mathbf{u, u})} = \\ = \sup\limits_{\|\mathbf{u}\| \ne 0} \frac{({\mathbf{Bu, Bu}})}{(\mathbf{u,u})}\sup\limits_{\|\mathbf{u}\| \ne 0}\frac{({\mathbf{B}}^{- 1}\mathbf{u},{\mathbf{B}}^{- 1}\mathbf{u})}{(\mathbf{u, u})} = {\|{\mathbf{B}}\|}_3^2 {\left\|{{\mathbf{B}}^{- 1}}\right\|}_3^2, $$$

откуда

$$\mu (\mathbf{A}) = \mu ^2({\mathbf{B}}).$$

Показать, что норма матрицы

$${\|\mathbf{A}\|}_2 = \max\limits_{1 \le j \le n}\sum\limits_{j= 1}^n{|{a_{ij}}|}$$

согласована с нормой вектора

$${\|\mathbf{u}\|}_2 = \sum\limits_{i = 1}^n{|{u_i}|}.$$

Решение.

$${\|\mathbf{Au}\|}_2 = \sum\limits_{i = 1}^n{\left|{\sum\limits_{j= 1}^n{a_{ij}u_j}}\right|} \le \sum\limits_{i = 1}^n{\sum\limits_{j = 1}^n{|{a_{ij}}|}} \cdot |{u_j}| \le \sum\limits_{j = 1}^n{|{u_j}|} \cdot \sum\limits_{i = 1}^n{|{a_{ij}}|} \le \\ \le (\max\limits_{1 \le i \le n}\sum\limits_j{|{a_{ij}}|}) \cdot {\|\mathbf{u}\|}_2 = {\|\mathbf{A}\|}_2 {\|\mathbf{u}\|}_2.$$

Положим

$$\max \sum\limits_{i = 1}^n{|{a_{ij}}|} = \sum\limits_{i = 1}^n{|{a_{ik}}|}.$$

Покажем, что существует вектор $$\mathbf{v}$$, для которого достигается равенство. В качестве такового можно взять вектор $$\mathbf{v}$$ с компонентами $${v}_i = 0, i \ne k, {v}_k = 1.$$

Таким образом, норма матрицы $${\|\mathbf{A}\|}_2 = \max\limits_{1 \le j \le n}\sum {|{u_{ij}}|}$$ согласована с нормой вектора $${\|\mathbf{u}\|}_2 = \sum\limits_{i = 1}^n{|{u_i}|}.$$

Дана жорданова клетка порядка n

$$\mathbf{A}= \left( \begin{array}{ccccccc} 1 {d} 0 \ldots \ldots 0 0 \\ 0 1 {d} \ldots \ldots 0 0 \\ \ldots \ldots \ldots \ldots \ldots \ldots \ldots \\ 0 0 \ldots 1 {d} \ldots 0 \\ \ldots \ldots \ldots \ldots \ldots \ldots \ldots \\ 0 0 \ldots \ldots \ldots 1 {d}\\ 0 0 \ldots \ldots \ldots \ldots 1 \\ \end{array} \right).$$

Найти $$\mu (\mathbf{A})$$ и оценить возмущение в компоненте $$u_1$$ решения системы $$\mathbf{Au} = \mathbf{f}$$, если компонент fn вектора f возмущен на величину $$\varepsilon.$$

Решение. Из $$\mathbf{Au}= \mathbf{f}$$ следует, что $$\mathbf{u}=\mathbf{A}^{- 1}\mathbf{f}.$$

С помощью обратной подстановки un = 1, un - 1 = ... находим компоненты матрицы

$$\mathbf{A}^{- 1}= \left( \begin{array}{ccccccc} 1 {- d} {d^2} \ldots \ldots {( - d)^{n - 2}} {( - d)^{n - 1}}\\ 0 1 {- d} \ldots \ldots {( - d)^{n - 3}} {( - d)^{n - 2}}\\ 0 0 1 \ldots \ldots {( - d)^{n - 4}} {( - d)^{n - 3}}\\ \ldots \ldots \ldots \ldots \ldots \ldots \ldots \\ \ldots \ldots \ldots \ldots \ldots \ldots \ldots \\ 0 0 0 \ldots \ldots 1 {- d}\\ 0 0 0 \ldots \ldots 0 1 \\ \end{array} \right).$$

В этом случае $$$ {\|\mathbf{A}\|}_1 = 1 + |{d}|, {\left\|{\mathbf{A}^{- 1}}\right\|}_1 = 1 + |{d}|+{d}^2 + \ldots + {|{d}|}^{n - 1} = \frac{{|{d}|}^n - 1}{|{d}| - 1} $;$$

Видно, что при | d | > 1 матрица $$\mathbf{A}$$ плохо обусловлена, при | d | < 1 - хорошо. При n = 20 и d = 5 имеем $$\mu (\mathbf{A}) \approx 10^{14}.$$

Компонент $$$ \tilde{u}_1 $$$ решения возмущенной системы $$\tilde{\mathbf{u}} = {\mathbf{A}}^{- 1}{\mathbf{\tilde{f}}}$$ будет

$$$ \tilde{\mathbf{u}}_1 = f_1 - df_2 + d^2f_3 + \ldots + {(- d)}^{n - 2}f_{n - 1} + {( - d)}^{n - 1}(f_n + \varepsilon ) = u_1 + {(- d)}^{n - 1}\varepsilon, $$$

где u1 — компонент решения невозмущенной системы $$\mathbf{Au} = \mathbf{f}$$

Отсюда видно, что при | d | > 1 возмущение в n компоненте вектора $$\mathbf{f}$$ увеличивается в компоненте u1 вектора $$\mathbf{u}$$ в | d |n - 1 раз, а при | d | < 1 — в | d |n - 1 раз убывает.

Пусть в системе линейных уравнений

u1  + 0,99u2  = f1, 
0,99u1  + u2  = f2

вектор $$\mathbf{f} = {(f_1, f_2)}^{T}$$ получает приращение $$\mathbf{\Delta f} = {(\delta f_1, \delta f_2)}^{T}$$, а решение получает приращение $$\Delta\mathbf{u} = {(\delta u_1, \delta u_2)}^{T}$$, Найти наименьшее число $$\mu,$$ при котором независимо от $$\mathbf{f}$$ и $$\Delta \mathbf{f}$$ выполняется оценка$$\frac{\|{\Delta\mathbf{u}}\|}{\|\mathbf{u}\|} \le \mu \frac{\|{\Delta \mathbf{f}}\|}{\|\mathbf{f}\|}.$$

Решить задачу, используя нормы матриц $${\|\cdot\|}_1, {\|\cdot\|}_2 , {\|\cdot\|}_3.$$

Решение. Для возмущенной задачи $$\mathbf{A}(\mathbf{u} + \Delta\mathbf{u}) = \mathbf{f} + \Delta\mathbf{f}$$, из линейности системы следует $$\mathbf{A}\Delta\mathbf{u} = \Delta\mathbf{f}.$$ Для возмущения решения выполняется равенство $$\Delta\mathbf{u} = {\mathbf{A}}^{- 1}\Delta\mathbf{f}.$$ Тогда $$$ \|{\Delta\mathbf{u}}\| \le \|{\mathbf{A}^{- 1}}\| \cdot \|\mathbf{A}\| \cdot \frac{\|{\Delta\mathbf{f}}\|}{\|\mathbf{f}\|} \cdot \frac{\|\mathbf{f}\|}{\|\mathbf{A}\|} $.$$

Отсюда сразу следует $$$ \|{\Delta\mathbf{u}}\| \le \left\|{{\mathbf{A}}^{- 1}}\right\| \cdot \|\mathbf{A}\| \frac{\|{\Delta\mathbf{f}}\|}{\|\mathbf{f}\|}\|\mathbf{u}\| $$$, так как $$$ \frac{\|\mathbf{f}\|}{\|\mathbf{A}\|} \le \|{\Delta\mathbf{u}}\| $.$$

Тогда искомая оценка будет $$$ \frac{\|{\Delta\mathbf{u}}\|}{\|\mathbf{u}\|} \le \left\|{{\mathbf{A}}^{- 1}}\| \cdot \|\mathbf{A}\| \frac{\|{\Delta\mathbf{f}}\|}{\|\mathbf{f}\|}.$$ Обозначим $$\mu = \|{{\mathbf{A}}^{- 1}}\| \cdot \|\mathbf{A}\| $.$$ В этом случае наименьшим числом, при котором выполняется оценка $$$ \frac{\|{\Delta\mathbf{u}}\|}{\|\mathbf{u}\|} \le \mu \frac{\|{\Delta\mathbf{f}}\|}{\|\mathbf{f}\|} $$$, является $$\mu = \|{{\mathbf{A}}^{- 1}}\| \cdot \|\mathbf{A}\|.$$ Это — число обусловленности системы уравнений. Численное решение в соответствующих нормах получается легко.

При заданном фиксированном $$\mathbf{f}$$ найти наименьшее число $$\nu,$$ при котором независимо от $$\mathbf{\Delta f}$$ выполняется оценка

$$$ \frac{\|{\Delta\mathbf{u}}\|}{\|\mathbf{u}\|} \le \nu (f) \frac{\|{\Delta\mathbf{f}}\|}{\|\mathbf{f}\|} $.$$

Найти такую правую часть системы $$\mathbf{f}$$, которой соответствует наименьшее $$\nu,$$ а также само это значение при использовании третьей нормы матрицы.

Решение.

По условию задачи $$$ \nu \ge \frac{\|\mathbf{f}\|}{\|\mathbf{u}\|} \cdot \frac{\|{\Delta\mathbf{u}}\|}{\|{\Delta\mathbf{f}}\|} $.$$ Рассмотрим, какие значения может принимать это число. Точная нижняя грань для такой оценки, очевидно, $$$ \inf\limits_{\|{\Delta\mathbf{f}}\| \ne 0}\nu = \frac{\|\mathbf{f}\|}{\|\mathbf{u}\|} \cdot \sup\limits_{\|{\Delta\mathbf{f}}\| \ne 0} \frac{\left\|{{\mathbf{A}}^{- 1} \cdot \Delta\mathbf{f}}\right\|}{\|{\Delta\mathbf{f}}\|} $.$$ Так как надо найти оценку, не зависящую от начального возмущения (при решении конкретной задачи оно, очевидно, неизвестно), получим

$$$ \nu = \inf\limits_{\|{\Delta\mathbf{f}}\| \ne 0}\nu = {\left\|{{\mathbf{A}}^{- 1}}\right\|} \cdot \frac{\|\mathbf{f}\|}{\|\mathbf{u}\|}. $$$

Для точной нижней грани выполнено $$$ \inf\limits_\mathbf{f}\nu (\mathbf{f}) = \left\|{{\mathbf{A}}^{- 1}}\right\| \times \\ \times {\left({\sup\frac{\|{{\mathbf{A}}^{- 1}\mathbf{f}}\|}{\|\mathbf{f}\|}}\right)}^{- 1} = 1 $.$$ Можно оценить и точную верхнюю грань:

$$$ \sup\limits_\mathbf{f}\nu = \left\|{{\mathbf{A}}^{- 1}}\right\| \sup\limits_{\mathbf{u}}\frac{\|\mathbf{Au}\|}{\|\mathbf{u}\|} = \left\|{{\mathbf{A}}^{- 1}}\right\| \cdot \|\mathbf{A}\| = \mu $$$

Таким образом, $$1 \le \nu \le \mu.$$

Ответим на вопрос, при каких $$\mathbf{f}$$ достигается $$\sup\limits_{\mathbf{f}}\nu $$ и $$\inf\limits_{\mathbf{f}}\nu.$$ Для этого используем разложение вектора правой части системы по базису из собственных векторов матрицы $${\mathbf{A}}^2$$ (без ограничения общности полагаем, что такой базис существует). В этом базисе$$\frac{\|\mathbf{Au}\|}{\|\mathbf{u}\|} = \sqrt{\frac{({\mathbf{A}}^*\mathbf{Au, u})}{(\mathbf{u, u})}} - \sqrt{\frac{(\sum{\lambda_i\xi_i\omega_i, \xi_i\omega_i})}{(\xi_i\omega_i, \xi_i\omega_i )}} = \sqrt{\frac{\sum{\lambda_i\xi_i^2}}{\sum{\xi_i^2}}},$$ $$$ \sqrt{\lambda_{\min}({\mathbf{A}}^*\mathbf{A})} \le \sqrt{\frac{\sum {\lambda_i\xi_i^2}}{\sum{\xi_i^2}}} = \frac{\|\mathbf{Au}\|}{\|\mathbf{u}\|} \le \sqrt{\lambda_{\max}({\mathbf{A}}^*\mathbf{A})} $.$$

Для самосопряженной положительной матрицы $${\mathbf{A}}^* = \mathbf{A} > 0$$ получаем

$$$ |{\lambda_{\min}(\mathbf{A})}| \le \sqrt{\frac{\sum\limits_i {\lambda_i \xi_i^2}}{\sum\limits_i{\xi_i^2}}} = \frac{\|\mathbf{Au}\|}{\|\mathbf{u}\|} \le |{\lambda_{\max}(\mathbf{A})}|. $$$

В то же время

$$$ \frac{\|{{\mathbf{A\omega}}_i}\|}{\|{\omega_i}\|} = \sqrt{\frac{(\mathbf{A}*{\mathbf{A\omega}}_i, \omega_i )}{(\omega_i, \omega_i )}} = \sqrt{\frac{\lambda_i (\omega_i, \omega_i )}{(\omega_i ,\omega_i )}} = \sqrt{\lambda_i (\mathbf{A}*\mathbf{A})}, $$$

для самосопряженной положительной матрицы $$$ \frac{\|{{\mathbf{A\omega}}_i}\|}{\|{\omega_i}\|} = |{\lambda_i (\mathbf{A})}| $.$$

$$$ \sup\limits_\omega\frac{\|{{\mathbf{A\omega}}}\|}{\|\omega\|} = \sqrt{\lambda_{\max}({\mathbf{A}}^*\mathbf{A})};\quad \inf\limits_\omega \frac{\|{{\mathbf{A\omega}}}\|}{\|\omega\|} = \sqrt{\lambda_{\min}({\mathbf{A}}^*\mathbf{A})}. $$$

В случае $${\mathbf{A}}^* = \mathbf{A} > 0$$ получаем $$\sup\limits_\omega \frac{\|{{\mathbf{A\omega}}}\|}{\|\omega\|} = |{\lambda_{\max}(\mathbf{A})}| \inf\limits_\omega \frac{\|{{\mathbf{A\omega}}}\|}{\|\omega\|} = |{\lambda_{\min}(\mathbf{A})}| $.$$

Таким образом,

$$$ \sup\limits_\mathbf{u} \frac{\|\mathbf{Au}\|}{\|\mathbf{u}\|} = \sqrt{\lambda_{\max}({\mathbf{A}}^*\mathbf{A});} \sup\limits_\omega\frac{\|{{\mathbf{A\omega}}}\|}{\|\omega\|} = \sqrt{\lambda_{\max}({\mathbf{A}}^*\mathbf{A})}, $$$

$$$ \inf\limits_\mathbf{u}\frac{\|\mathbf{Au}\|}{\|\mathbf{u}\|} = \sqrt{\lambda_{\min}({\mathbf{A}}^*\mathbf{A});} \inf\limits_\omega \frac{\|{{\mathbf{A\omega}}}\|}{\|\omega\|} = \sqrt{\lambda_{\min}({\mathbf{A}}^*\mathbf{A})}. $$$

Выписать формулы итерационных методов Якоби, Зейделя, верхней релаксации для СЛАУ

$$\mathbf{Ax}= \mathbf{f}, \mathbf{A}= \left( \begin{array}{cc} 2 1 \\ 1 2 \\ \end{array}\right), \mathbf{x} = {(1,-1)}^{T}, \mathbf{f} = {(1,-1)}^{T}$$

или

2u + v = 1, 
u + 2v = - 1.

Оценить количество итераций для метода Якоби.

(Решение системы: u = 1; v = -1 ).

Решение.

Итерационные методы Якоби, Зейделя, релаксации соответственно записываются

$$\left\{ \begin{array}{ccc} u_{k + 1} = - \frac{1}{2}v_k + \frac{1}{2}, \\ v_{k + 1} = - \frac{1}{2}u_k - \frac{1}{2}, \\ \end{array} \right.$$ $$(u_0, v_0 ) = (u_0^0, v_0^0 )$$,

или $${\mathbf{x}}_{k+ 1} = {\mathbf{Bx}}_k + \mathbf{f}$$,

$${\mathbf{B}}= \left( \begin{array}{cc} 0 {- 0,5}\\ {- 0,5} 0 \\ \end{array} \right),$$

$$\begin{gather*} u_{k + 1} = - \frac{1}{2}v_k + \frac{1}{2}, \\ v_{k + 1} = - \frac{1}{2}u_{k + 1} - \frac{1}{2}; \\ u_{k + 1} = (1 - \tau )u_k + \frac{\tau}{2}(1 - v_k), \\ v_{k + 1} = (1 - \tau )v_k - \frac{\tau}{2}(1 + u_k). \end{gather*} $$

Оценка количества итераций проводится по формуле

$$$ k \approx \ln \frac{\varepsilon }{\varepsilon_0 }/\ln \|\mathbf{B}\| = \ln 10^{- 3}/\ln \frac{1}{2}. $$$

Представить графическую интерпретацию итерационного метода Якоби для СЛАУ

$$\begin{gather*} a_{11}u + a_{12}v = f_1, \\ a_{21}u + a_{22}v= f_2, \\ a_{11} \ne 0, a_{22} \ne 0 \end{gather*}$$

Решение. Итерационный процесс Якоби записывается как

$$\begin{gather*}u_{k+ 1} = - \frac{a_{12}}{a_{11}}v_k + \frac{f_1}{a_{11}}, \\ v_{k+ 1} = - \frac{a_{21}}{a_{22}}u_k + \frac{f_2}{a_{22}}. \end{gather*}$$(рис 2.1)

Первое уравнение соответствует прямой 1, второе — прямой 2. Вычисление u1 соответствует проведению отрезка, параллельного оси 0u и (при v = v0 ) до пересечения с прямой 1 ; точка пересечения даст первое приближение u1. Вычислению v1 соответствует проведение из точки A1 прямой, параллельной оси 0v до пересечения с прямой 2 и т.д. до сходимости итераций к точке пересечения прямых 1 и 2 (Ak) с заданной точностью.

При каких a,b сходится метод простой итерации $${\mathbf{u}}_{k+ 1} = {\mathbf{Bu}}_k + \mathbf{f}$$, где

$${\mathbf{B}} = \left( \begin{array}{ccc} {a} {b} 0 \\ {b} {a} {b}\\ 0 {b} {a} \\ \end{array} \right).$$

Решение. Для того, чтобы метод простой итерации сходился к решению соответствующей СЛАУ, необходимо и достаточно, чтобы все собственные значения матрицы $$\mathbf{B}$$ по модулю были меньше единицы: $$|\lambda _{i}| < 1.$$ Решаем характеристическое уравнение

$$\begin{gather*} det (\mathbf{B}- \lambda \mathbf{E}) = {\left| \begin{array}{ccc} {a - \lambda } {b} 0 \\ {b} {{a}-{\lambda}} {b}\\ 0 {b} {{a}-{\lambda}}\\ \end{array} \right|} = \\ = (a - \lambda ) {\left| \begin{array}{cc} {{a}-{\lambda}} {b}\\ {b} {{a}-{\lambda}}\\ \end{array} \right|} - b{\left| \begin{array}{cc} {b} {a}\\ 0 {{a}-{\lambda}}\\ \end{array} \right|} = \\ = (a - \lambda )\left[{{(a - \lambda )}^2 - b^2}\right] - b^2 (a - \lambda ) = \\ = (a - \lambda )(a - \lambda - \sqrt{2b})(a - \lambda + \sqrt{2b}) = 0, \end{gather*}$$

откуда получим условие сходимости итерационного метода

$$|{a}| < 1, |{a \pm \sqrt{2b}}| < 1.$$

Найти условие сходимости итерационных методов Якоби и Зейделя для СЛАУ $$\mathbf{Au} = \mathbf{f}$$ с матрицей

$$\mathbf{A}$$ вида$$\mathbf{A} = \left( \begin{array}{ccc} {a} {b} 0 \\ {b} {a} {b}\\ 0 {b} {a}\\ \end{array} \right).$$

Решение. Для метода Якоби

$${\mathbf{u}}_{k + 1} = {\mathbf{Bu}}_k + \mathbf{f},\quad \mbox{где}\quad {\mathbf{B}} = -{\mathbf{D}}^{- 1}({\mathbf{L}} + {\mathbf{U}}).$$

Имеет место уравнение:

$${\mathbf{B\omega}} = \lambda \omega$$, где $$\lambda$$ и $$\mathbf{\omega}$$ — собственное число и собственный вектор, соответственно. В таком случае $$- {\mathbf{D}}({\mathbf{L}} + {\mathbf{U}})\omega = \lambda\omega$$, или: $$({\mathbf{L}} + {\mathbf{U}} + \lambda{\mathbf{D}}) \omega = 0$$, откуда (предполагаем наличие нетривиальных решений у последней СЛАУ ):

$$\det ({\mathbf{L}} + {\mathbf{U}} + \lambda{\mathbf{D}}) = 0.$$

Решим это уравнение:

$$\det {\left| \begin{array}{ccc} \lambda a b 0\\ b \lambda a b\\ 0 b \lambda a \end{array} \right|} = 0,$$

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

$$$ |{\frac{b}{a}}| < 2^{- \frac{1}{2}}. $$$

Для метода Зейделя имеем

$${\mathbf{B}} = {- (\mathbf{D} + \mathbf{L})}^{- 1}{\mathbf{U}},\quad {\mathbf{B\omega}} = \lambda\omega.$$

В таком случае

$${- (\mathbf{D} + \mathbf{L})}^{- 1}{\mathbf{U\omega}} = \lambda\omega$$,

откуда следует уравнение

$$\det (\lambda{\mathbf{L}} + \lambda {\mathbf{D}} + {\mathbf{U}}) = 0$$,

Вычислив детерминант, придем к алгебраическому уравнению

$$\det {\left( \begin{array}{ccc} \lambda a b 0 \\ \lambda b \lambda a b \\ 0 \lambda b \lambda a \end{array} \right) } = a{\lambda}^2(a^2\lambda - 2b^2) = 0.$$

В таком случае, поскольку $$$ \lambda_{1,2}= 0, \lambda_3 = 2\frac{b^2}{a^2} $$$, получим условие сходимости метода Зейделя: $$$ \frac{b}{a} < 2^{- \frac{1}{2}}$.$$ Видно, что в данном случае условия сходимости для обоих методов совпадают.