Якоби, Карл Густав Якоб
- 1 year ago
- 0
- 0
Метод Якоби — разновидность метода простой итерации для численного решения системы линейных алгебраических уравнений . Назван в честь Карла Густава Якоби .
Пусть требуется численно решить систему линейных уравнений:
Предполагается, что (иначе метод Якоби неприменим). Выразим через первое уравнение, — через второе и т. д. :
В методе Якоби последовательность приближений строится следующим образом. Выбирается первое приближение , формула для остальных приближений имеет вид
В матричной форме имеет следующий вид. Пусть СЛАУ в матричной форме записано как
Представим матрицу в виде , где означает диагональную матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы ; тогда как матрицы и содержат верхнюю и нижнюю треугольные части , на главной диагонали которых нули. Тогда итерационную формулу можно записать как
Приведем достаточное условие сходимости метода.
Теорема
.
Пусть . Тогда при любом выборе начального приближения :
|
Условие окончания итерационного процесса при достижении точности имеет вид:
Для достаточно хорошо обусловленной матрицы (при ) оно выполняется при
Данный критерий не требует вычисления нормы матрицы и часто используется на практике. При этом точное условие окончания итерационного процесса имеет вид
и требует дополнительного умножения матрицы на вектор на каждой итерации, что примерно в два раза увеличивает вычислительную сложность алгоритма.
В отличие от метода Гаусса-Зейделя мы не можем заменять на в процессе итерационной процедуры, так как эти значения понадобятся для остальных вычислений. Это наиболее значимое различие между методом Якоби и методом Гаусса-Зейделя решения СЛАУ . Таким образом на каждой итерации придётся хранить оба вектора приближений: старый и новый.
Ниже приведён алгоритм реализации на 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