№11 Решение интегрального уравнения Фредгольма в случае вырожденного ядра

eps = 1e-6
V = 8.0

def f(x):
    return V * (4 / 3 * x + 1 / 4 * x ** 2 + 1 / 5 * x ** 3)

def check_y(x):
    return V * x

def gauss(sole):
    n = len(sole)
    m = len(sole[0])
    x_coords = list(range(n))

    for j in range(m - 1):
        mx = j
        for i in range(j, n):
            if abs(sole[i][j]) > abs(sole[mx][j]):
                mx = i

        if abs(sole[mx][j]) < eps:
            continue

        sole[j], sole[mx] = sole[mx], sole[j]

        mx = j
        for i in range(j, n):
            if abs(sole[j][i]) > abs(sole[j][mx]):
                mx = i

        x_coords[j], x_coords[mx] = x_coords[mx], x_coords[j]
        for i in range(n):
            sole[i][j], sole[i][mx] = sole[i][mx], sole[i][j]

        if abs(sole[j][mx]) < eps:
            raise ValueError("Division by zero!")

        for i in range(j + 1, n):
            d = sole[i][j] / sole[j][j]
            for k in range(m):
                sole[i][k] -= sole[j][k] * d

    for j in range(n - 1, -1, -1):
        for i in range(j - 1, -1, -1):
            d = sole[i][j] / sole[j][j]
            sole[i][j] -= sole[j][j] * d
            sole[i][m - 1] -= sole[j][m - 1] * d

    s = [0] * n
    for i in range(n):
        s[x_coords[i]] = sole[i][m - 1] / sole[i][i]

    return s

def gauss_with_b(A, B):
    for i in range(len(A)):
        A[i].append(B[i])
    return gauss(A)

def print_results(vx, vy):
    k = 5
    for i in range(0, len(vx), k):
        m = min(i + k, len(vx))
        print("\t".join(f"{vx[j]:.5f}" for j in range(i, m)))
        print("\t".join(f"{vy[j]:.5f}" for j in range(i, m)))
        print("\t".join(f"{check_y(vx[j]):.5f}" for j in range(i, m)))
        print("\t".join(f"{vy[j] - check_y(vx[j]):.5f}" for j in range(i, m)))
        print()

def solve():
    n = 3
    A = [[1 / (i + j + 3) for j in range(n)] for i in range(n)]
    for i in range(n):
        A[i][i] += 1

    B = [4 * V / 3 / (i + 3) + V / 4 / (i + 4) + V / 5 / (i + 5) for i in range(n)]
    vq = gauss_with_b(A, B)

    m = 3
    l, r = 0, 1
    h = (r - l) / (m - 1)
    vx = [l + i * h for i in range(m)]

    vy = [f(x) - sum(vq[k] * x ** (k + 1) for k in range(n)) for x in vx]
    print_results(vx, vy)


if __name__ == "__main__":
    solve()