Русская Википедия:Метод Якоби
Шаблон:Об Шаблон:Другие значения термина Метод Якоби — разновидность метода простой итерации для численного решения системы линейных алгебраических уравнений. Назван в честь Карла Густава Якоби.
Описание метода
Пусть требуется численно решить систему линейных уравнений:
<math>\left\{ \begin{array}{rcl} a_{11}x_1 + \ldots + a_{1n}x_n& = & b_{1} \\ \ \ldots 14:14, 27 августа 2023 (+04)14:14, 27 августа 2023 (+04)14:14, 27 августа 2023 (+04)14:14, 27 августа 2023 (+04)EducationBot (обсуждение) \\ 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:
- <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
Условие остановки
Условие окончания итерационного процесса при достижении точности <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++
#include <math.h>
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
Примечания
Литература
См. также
Шаблон:Rq Шаблон:Методы решения СЛАУ
