Метод Якоби

Шаблон:Об Шаблон:Другие значения термина Метод Якоби — разновидность метода простой итерации для численного решения системы линейных алгебраических уравнений. Назван в честь Карла Густава Якоби.

Описание метода

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

$\{\begin{array}{ccc}{a}_{11}{x}_1+\dots +a_{1n}{x}_n& =& b_1\\ \dots \\ a_{n1}{x}_1+\dots +a_{nn}{x}_n& =& b_n\end{array}$

Предполагается, что $a_{ii}\ne 0,i\in \{1,\dots ,n\}$ (иначе метод Якоби неприменим). Выразим $x_1$ через первое уравнение, $x_2$ — через второе и т. д.Шаблон:Sfn:

$\{\begin{array}{cc}{x}_1=& \frac{1}{a_{11}}(b_1-a_{12}{x}_2-\dots -a_{1n}{x}_n)\\ \dots \\ x_n=& \frac{1}{a_{nn}}(b_n-a_{n1}{x}_1-\dots -a_{n(n-1)}{x}_{n-1})\end{array}$

В методе Якоби последовательность приближений $x^{(k)}$ строится следующим образом. Выбирается первое приближение $x^{(0)}$, формула для остальных приближений имеет вид

$\begin{array}{cc}{x}_1^{(k+1)}=& \frac{1}{a_{11}}(b_1-a_{12}{x}_2^{(k)}-\dots -a_{1n}{x}_n^{(k)})\\ \dots \\ x_n^{(k+1)}=& \frac{1}{a_{nn}}(b_n-a_{n1}{x}_1^{(k)}-\dots -a_{n(n-1)}{x}_{n-1}^{(k)})\end{array}$.

В матричной форме имеет следующий вид. Пусть СЛАУ в матричной форме записано как

$Ax=b$

Представим матрицу $A$ в виде $A=L+D+U$, где $D$ означает диагональную матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы $A$; тогда как матрицы $U$ и $L$ содержат верхнюю и нижнюю треугольные части $A$, на главной диагонали которых нули. Тогда итерационную формулу можно записать как

$x^{(k+1)}=D^{-1}(b-(L+U)x^{(k)})$

Сходимость

Приведем достаточное условие сходимости метода. Шаблон:Message box

Условие остановки

Условие окончания итерационного процесса при достижении точности $\epsilon$ имеет вид:

$\parallel x^{(k+1)}-x^{(k)}\parallel \le \frac{1-q}{q}\epsilon .$

Для достаточно хорошо обусловленной матрицы $B$ (при $\parallel B\parallel =q<1/2$) оно выполняется при

$\parallel x^{(k+1)}-x^{(k)}\parallel \le \epsilon .$

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

$\parallel A{x}^{(k)}-b\parallel \le \epsilon$

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

Сравнение с другими методами

В отличие от метода Гаусса-Зейделя мы не можем заменять $x_i^{(k)}$ на $x_i^{(k+1)}$ в процессе итерационной процедуры, так как эти значения понадобятся для остальных вычислений. Это наиболее значимое различие между методом Якоби и методом Гаусса-Зейделя решения СЛАУ. Таким образом на каждой итерации придётся хранить оба вектора приближений: старый и новый.

Реализация

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

#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;
}

Ниже приведен алгоритм реализации на 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

Литература

• Шаблон:Книга

См. также

• Метод
• Метод Гаусса — Зейделя