Решение системы уравнений методами Якоби, Зейделя, Гаусса

import numpy as np
 
def jacobi(A, b, x0, eps=1e-5, max_iter=100):
    D = np.diag(A)
    R = A - np.diagflat(D)
    x = x0
    for i in range(max_iter):
        x_new = (b - np.dot(R, x)) / D
        if np.allclose(x, x_new, rtol=eps):
            return x_new
        x = x_new
    return x
 
# Пример использования
A = np.array([[10, -1, 2, 0],
              [-1, 11, -1, 3],
              [2, -1, 10, -1],
              [0, 3, -1, 8]])
b = np.array([6, 25, -11, 15])
x0 = np.zeros(4)
 
x = jacobi(A, b, x0)
print(x)

import numpy as np
 
def gauss_seidel(A, b, x0, eps=1e-5, max_iter=100):
    D = np.diag(A)
    L = np.tril(A, k=-1)
    U = np.triu(A, k=1)
    x = x0
    for i in range(max_iter):
        x_new = (b - np.dot(L, x) - np.dot(U, x)) / D
        if np.allclose(x, x_new, rtol=eps):
            return x_new
        x = x_new
    return x
 
# Пример использования
A = np.array([[10, -1, 2, 0],
              [-1, 11, -1, 3],
              [2, -1, 10, -1],
              [0, 3, -1, 8]])
b = np.array([6, 25, -11, 15])
x0 = np.zeros(4)
 
x = gauss_seidel(A, b, x0)
print(x)

import numpy as np
 
def gauss(A, b):
    n = len(A)
 
    for i in range(n):
        # Находим максимальный элемент в столбце
        max_elem = abs(A[i][i])
        max_row = i
        for k in range(i + 1, n):
            if abs(A[k][i]) > max_elem:
                max_elem = abs(A[k][i])
                max_row = k
 
        # Переставляем строки матрицы, если требуется
        if i != max_row:
            A[i], A[max_row] = A[max_row], A[i]
            b[i], b[max_row] = b[max_row], b[i]
 
        # Обнуляем элементы ниже главной диагонали
        for k in range(i + 1, n):
            c = -A[k][i] / A[i][i]
            for j in range(i, n):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]
            b[k] += c * b[i]
 
    # Находим решение
    x = [0 for i in range(n)]
    for i in range(n - 1, -1, -1):
        x[i] = b[i] / A[i][i]
        for k in range(i - 1, -1, -1):
            b[k] -= A[k][i] * x[i]
    return x
 
# Пример использования
A = np.array([[10, -1, 2, 0],
              [-1, 11, -1, 3],
              [2, -1, 10, -1],
              [0, 3, -1, 8]])
b = np.array([6, 25, -11, 15])
 
x = gauss(A, b)
print(x)