Просмотр кода страницы Метод Якоби

{{Об|методе решения систем линейных алгебраических уравнений|итерационном алгоритме нахождения собственных значений|метод Якоби для собственных значений}}
{{другие значения термина|Якоби|Якоби}}
'''Метод Якоби''' — разновидность [[простой итерации метод|метода простой итерации]] для численного решения [[Система линейных алгебраических уравнений|системы линейных алгебраических уравнений]]. Назван в честь [[Якоби, Карл Густав Якоб|Карла Густава Якоби]].

== Описание метода ==
Пусть требуется численно решить систему линейных уравнений:

<math>\left\{
\begin{array}{rcl}
a_{11}x_1 + \ldots + a_{1n}x_n& = & b_{1} \\
\ \ldots ~~~~~~~~~~~~~~~~~~~~~~~ \\
a_{n1}x_1 + \ldots + a_{nn}x_n & = & b_{n} 
\end{array} \right.</math>

Предполагается, что <math>a_{ii} \neq 0, \ i \in \{1,\ldots,n\}</math> (иначе метод Якоби неприменим). Выразим <math>x_1</math> через первое уравнение, <math>x_2</math> — через второе и т. д.{{sfn|Березин, Жидков|1959|с=57}}:
:<math>\left\{
\begin{array}{rcl}
x_1 = & \dfrac{1}{a_{11}}\left(b_{1} - a_{12}x_2 - \ldots - a_{1n}x_n\right) \\
\ldots  \\
x_n = & \dfrac{1}{a_{nn}}\left(b_{n} - a_{n1}x_1 - \ldots - a_{n(n-1)}x_{n-1}\right)
\end{array} \right.</math>

В методе Якоби последовательность приближений <math>x^{(k)}</math> строится следующим образом. Выбирается первое приближение <math>x^{(0)}</math>, формула для остальных приближений имеет вид
:<math>
\begin{array}{rcl}
x_1^{(k+1)} = & \dfrac{1}{a_{11}}\left(b_{1} - a_{12}x_2^{(k)} - \ldots - a_{1n}x_n^{(k)}\right) \\
\ldots  \\
x_n^{(k+1)} = & \dfrac{1}{a_{nn}}\left(b_{n} - a_{n1}x_1^{(k)} - \ldots - a_{n(n-1)}x_{n-1}^{(k)}\right)
\end{array}</math>.

В матричной форме имеет следующий вид. Пусть СЛАУ в матричной форме записано как
:<math>Ax=b</math>
Представим матрицу <math>A</math> в виде <math>A=L+D+U</math>, где <math>D</math> означает диагональную матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы <math>A</math>; тогда как матрицы <math>U</math> и <math>L</math> содержат верхнюю и нижнюю треугольные части <math>A</math>, на главной диагонали которых нули. Тогда итерационную формулу можно записать как
:<math>x^{(k+1)} = D^{-1}(b-(L+U)x^{(k)})</math>

== Сходимость ==
Приведем достаточное условие сходимости метода.
{{message box |backgroundcolor = white |image =Logo_arte.jpg |size = 50px| heading=[[Метод итерации|Теорема]].|message= Пусть <math>\| B \| < 1</math>. Тогда при любом выборе  начального приближения <math>\vec{x}^{(0)}</math>:
# метод сходится;
# скорость сходимости метода равна скорости сходимости [[Геометрическая прогрессия|геометрической прогрессии]] со знаменателем <math>q= \|B\|</math>;
# верна оценка погрешности: <math>\|\vec{x}^{(k)}-\vec{x}\| < q^k \, \|\vec{x}^{(0)}-\vec{x}\|</math>.}}

== Условие остановки ==
Условие окончания итерационного процесса при достижении точности <math>\varepsilon</math> имеет вид:
: <math>\parallel x^{(k+1)}-x^{(k)}\parallel \le \frac{1-q}{q}\varepsilon.</math>

Для достаточно хорошо обусловленной матрицы <math>B</math> (при <math>\parallel B\parallel = q < 1/2 </math>) оно выполняется при
: <math>\parallel x^{(k+1)}-x^{(k)}\parallel \le \varepsilon.</math>
Данный критерий не требует вычисления нормы матрицы и часто используется на практике. При этом точное условие окончания итерационного процесса имеет вид
: <math>\parallel Ax^{(k)}-b\parallel \le \varepsilon</math>
и требует дополнительного умножения матрицы на вектор на каждой итерации, что примерно в два раза увеличивает вычислительную сложность алгоритма.

== Сравнение с другими методами ==
В отличие от метода [[Метод Зейделя|Гаусса-Зейделя]] мы не можем заменять <math>x^{(k)}_i</math> на <math>x^{(k+1)}_i</math> в процессе итерационной процедуры, так как эти значения понадобятся для остальных вычислений. Это наиболее значимое различие между методом Якоби и методом Гаусса-Зейделя решения [[СЛАУ]]. Таким образом на каждой итерации придётся хранить оба вектора приближений: старый и новый.

== Реализация ==
Ниже приведён алгоритм реализации на C++

<source lang="cpp">
#include <cmath>
const double eps = 0.001; ///< желаемая точность 

// ..........................

/// N - размерность матрицы; A[N][N] - матрица коэффициентов, F[N] - столбец свободных членов,
/// X[N] - начальное приближение, ответ записывается также в X[N];
void Jacobi (int N, double** A, double* F, double* X)
{
	double* TempX = new double[N];
	double norm; // норма, определяемая как наибольшая разность компонент столбца иксов соседних итераций.

	do {
		for (int i = 0; i < N; i++) {
			TempX[i] = F[i];
			for (int g = 0; g < N; g++) {
				if (i != g)
					TempX[i] -= A[i][g] * X[g];
			}
			TempX[i] /= A[i][i];
		}
        norm = fabs(X[0] - TempX[0]);
		for (int h = 0; h < N; h++) {
			if (fabs(X[h] - TempX[h]) > norm)
				norm = fabs(X[h] - TempX[h]);
			X[h] = TempX[h];
		}
	} while (norm > eps);
	delete[] TempX;
}
</source>Ниже приведен алгоритм реализации на [[Python]]<syntaxhighlight lang="python">
from collections.abc import Sequence, MutableSequence


def Jacobi(
        A: Sequence[Sequence[float]],
        b: Sequence[float],
        eps: float = 0.001,
        x_init: MutableSequence[float] | None = None) -> list[float]:
    """
    метод Якоби для A*x = b (*)

    :param A: Матрица (*) слева

    :param b: известный вектор (*) справа

    :param x_init: начальное приближение

    :return: приблизительное значения х (*)
    """

    def sum(a: Sequence[float], x: Sequence[float], j: int) -> float:
        S: float = 0
        for i, (m, y) in enumerate(zip(a, x)):
            if i != j:
                S += m*y
        return S

    def norm(x: Sequence[float], y: Sequence[float]) -> float:
        return max((abs(x0-y0) for x0, y0 in zip(x, y)))

    if x_init is None:
        y = [a/A[i][i] for i, a in enumerate(b)]
    else:
        y = x.copy()

    x: list[float] = [-(sum(a, y, i) - b[i])/A[i][i]
                      for i, a in enumerate(A)]

    while norm(y, x) > eps:
        for i, elem in enumerate(x):
            y[i] = elem
        for i, a in enumerate(A):
            x[i] = -(sum(a, y, i) - b[i])/A[i][i]
    return x

</syntaxhighlight>

== Примечания ==
{{Примечания}}

== Литература ==
* {{книга |автор={{автор|Березин, Иван Семёнович|Березин, И. С.}}, {{автор|Жидков, Николай Петрович|Жидков Н. П.}} |заглавие=Методы вычислений |том=II |место=М. |издательство=[[Физматлит]] |год=1959 |язык=ru |ref=Березин, Жидков }}

== См. также ==
* [[Метод ]]
* [[Метод Зейделя|Метод Гаусса — Зейделя]]

{{rq|topic=math|sources}}
{{Методы решения СЛАУ}}


[[Категория:Методы решения СЛАУ]]