Участник:Sergey.protserov/Метод Якоби решения СЛАУ

Свойства и структура алгоритмов

Общее описание алгоритма

Метод Якоби -- одношаговый стационарный итерационный метод решения СЛАУ вида [math]\displaystyle{ Ay = f }[/math], где [math]\displaystyle{ A = \left( \begin{array}{ccc} a_{11} & \dots & a_{1m} \\ \vdots & \ddots & \vdots \\ a_{m1} & \dots & a_{mm} \\ \end{array} \right) }[/math], [math]\displaystyle{ f = \left( \begin{array}{c} f_{1} \\ \vdots \\ f_{m} \\ \end{array} \right) }[/math], [math]\displaystyle{ y = \left( \begin{array}{c} y_{1} \\ \vdots \\ y_{m} \\ \end{array} \right) }[/math], [math]\displaystyle{ \det A \ne 0 }[/math].

Каноническая форма одношагового стационарного итерационного метода имеет вид : [math]\displaystyle{ B\frac{y^{n+1} - y^{n}}{\tau} + Ay^{n} = f, \quad n = 0,\,1,\,\dots\,, }[/math]

где [math]\displaystyle{ B }[/math] — невырожденная матрица [math]\displaystyle{ m \times m }[/math], [math]\displaystyle{ \tau \in \mathbb{R} }[/math], [math]\displaystyle{ y^{0} }[/math] — заданное начальное приближение. Решение исходной СЛАУ находится приближённо посредством последовательных итераций. На [math]\displaystyle{ n }[/math]-ом шаге находится [math]\displaystyle{ y^{n+1} }[/math] — очередное приближение для искомого решения [math]\displaystyle{ y }[/math].

В методе Якоби [math]\displaystyle{ \tau = 1 }[/math], [math]\displaystyle{ B = D }[/math], где [math]\displaystyle{ D }[/math] — диагональная матрица, элементы которой совпадают с элементами, стоящими на главной диагонали матрицы [math]\displaystyle{ A }[/math].

Достаточным условием сходимости метода является свойство строгого диагонального преобладания у матрицы [math]\displaystyle{ A }[/math].

Математическое описание алгоритма

В обозначениях предыдущего пункта выражение для [math]\displaystyle{ y^{n+1} }[/math] через [math]\displaystyle{ y^{n} }[/math]: [math]\displaystyle{ y^{n+1} = D^{-1}\left(D-A\right)y^{n} + D^{-1}f }[/math].

В поэлементной записи:

[math]\displaystyle{ y^{n+1}_{i} = \frac{1}{a_{ii}}\left(f_{i} - \sum_{j=1,\,j \ne i}^{m}a_{ij}y^{n}_{j}\right),\quad i = 1,\,\dots,\,m }[/math].

В качестве условия окончания итерационного процесса можно использовать условие [math]\displaystyle{ \left\lVert y^{n+1} - y^{n}\right\rVert \le \varepsilon }[/math], где [math]\displaystyle{ \varepsilon }[/math] — заданная точность. Кроме того, можно ограничить максимальное число итераций, задав [math]\displaystyle{ n_{max} }[/math]. Для оценки ошибки можно использовать невязку [math]\displaystyle{ Ay^{n+1} - f }[/math].

Вычислительное ядро алгоритма

Основное время работы алгоритма приходится на последовательные вычисления векторов [math]\displaystyle{ y^{n+1} }[/math] по формуле, приведённой в предыдущем пункте, при уже вычисленных в начале работы алгоритма матрице [math]\displaystyle{ D^{-1}A }[/math] и векторе [math]\displaystyle{ D^{-1}f }[/math].

Макроструктура алгоритма

В описываемом алгоритме можно выделить следующие макрооперации:

1. вычисление [math]\displaystyle{ D^{-1} }[/math]
2. вычисление [math]\displaystyle{ D^{-1}A }[/math]
3. вычисление [math]\displaystyle{ D^{-1}f }[/math]
4. вычисление [math]\displaystyle{ y^{n+1} = \left(I - D^{-1}A\right)y^{n} + D^{-1}f }[/math]

Макрооперации 1-3 выполняются единожды, и в силу того, что матрица [math]\displaystyle{ D }[/math] — диагональная, занимают лишь незначительную часть времени работы алгоритма.

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

Кроме того, если в качестве критерия завершения работы алгоритма используется условие [math]\displaystyle{ \left\lVert y^{n+1} - y^{n}\right\rVert \le \varepsilon }[/math], требуется также вычислять указанную величину. В дальнейшем при описании алгоритма мы будем предполагать, что этот критерий не используется, а используется завершение работы по достижении максимального числа итераций.

Схема реализации последовательного алгоритма

1. составить диагональную матрицу [math]\displaystyle{ D }[/math]
2. вычислить [math]\displaystyle{ D^{-1} }[/math]
3. вычислить [math]\displaystyle{ D^{-1}A }[/math]
4. вычислить [math]\displaystyle{ D^{-1}f }[/math]
5. выполнять вычисления по формуле [math]\displaystyle{ y^{n+1} = \left(I - D^{-1}A\right)y^{n} + D^{-1}f }[/math], [math]\displaystyle{ n = 0,\,1,\,\dots,\,n_{max} }[/math]

При этом на [math]\displaystyle{ n }[/math]-ом шаге итераций необходимо хранить оба вектора [math]\displaystyle{ y^{n} }[/math], [math]\displaystyle{ y^{n+1} }[/math].

Последовательная сложность алгоритма

Шаг 2 предыдущего пункта требует выполнения [math]\displaystyle{ m }[/math] операций деления вещественных чисел, шаг 3 (с учётом того, что матрица [math]\displaystyle{ D^{-1} }[/math] — диагональная) требует выполнения [math]\displaystyle{ m^{2} }[/math] операций умножения вещественных чисел, аналогично шаг 4 требует [math]\displaystyle{ m }[/math] операций умножения, а каждая итерация шага 7 требует [math]\displaystyle{ m^{2} }[/math] умножений и [math]\displaystyle{ m^{2} + 2m }[/math] сложений/вычитаний.

Ресурс параллелизма алгоритма

Шаг 2 требует один ярус из [math]\displaystyle{ m }[/math] операций деления, шаг 3 требует BEGIN FIXME один ярус из [math]\displaystyle{ m^{2} }[/math] операций умножения, шаг 4 требует один ярус из [math]\displaystyle{ m }[/math] операций умножения, а каждая итерация шага 7 требует по [math]\displaystyle{ m }[/math] ярусов умножений и сложений (в каждом из ярусов — [math]\displaystyle{ m }[/math] операций) для выполнения умножения матрицы на вектор и ещё два яруса по [math]\displaystyle{ m }[/math] сложений/вычитаний END FIXME, причём итерации выполняются последовательно.

Входные и выходные данные алгоритма

Входные данные:

1. Вещественная [math]\displaystyle{ m \times m }[/math] матрица [math]\displaystyle{ A }[/math], вообще говоря, плотная
2. Вещественный [math]\displaystyle{ m }[/math]-мерный вектор правой части [math]\displaystyle{ f }[/math]
3. Вещественный [math]\displaystyle{ m }[/math]-мерный вектор начального приближения [math]\displaystyle{ y^{0} }[/math]
4. Максимальное число итераций алгоритма [math]\displaystyle{ n_{max} }[/math]

Выходные данные:

1. Вещественный [math]\displaystyle{ m }[/math]-мерный вектор приближённого решения [math]\displaystyle{ y^{n_{max}} }[/math]