№12 Квадратурный метод решения интегрального уравнения Фредгольма

import numpy as np

def f(x, v=8):
    return v * (4 / 3 * x + 1 / 4 * x ** 2 + 1 / 5 * x ** 3)

def gauss(A, B):
    n = len(A)
    A = np.array(A, dtype=float)
    B = np.array(B, dtype=float)
    for i in range(n):
        # Partial pivoting
        max_row = max(range(i, n), key=lambda r: abs(A[r, i]))
        if abs(A[max_row, i]) < 1e-6:
            raise ValueError("Division by zero detected in Gaussian elimination.")
        if max_row != i:
            A[[i, max_row]] = A[[max_row, i]]
            B[[i, max_row]] = B[[max_row, i]]

        # Eliminate column i
        for j in range(i + 1, n):
            factor = A[j, i] / A[i, i]
            A[j, i:] -= factor * A[i, i:]
            B[j] -= factor * B[i]

    # Back substitution
    x = np.zeros_like(B)
    for i in range(n - 1, -1, -1):
        x[i] = (B[i] - np.dot(A[i, i + 1:], x[i + 1:])) / A[i, i]
    return x

def main():
    v = 8
    n = 3

    l, r = 0, 1
    h = (r - l) / (n - 1)
    vx = np.linspace(l, r, n)

    A = np.zeros((n, n))
    F = np.zeros(n)

    for i in range(n):
        x_i = vx[i]
        for j in range(n):
            x_j = vx[j]
            A[i, j] = h * (x_j * x_i + x_j ** 2 * x_i ** 2 + x_j ** 3 * x_i ** 3)
            if i == j:
                A[i, j] += 1
        F[i] = f(x_i, v)

    y = gauss(A, F)

    res1 = np.zeros((4, n))
    res1[0] = vx

    for i in range(n):
        x_i = vx[i]
        res1[1, i] = f(x_i, v)
        for j in range(n):
            x_j = vx[j]
            res1[1, i] -= h * (x_j * x_i + x_j ** 2 * x_i ** 2 + x_j ** 3 * x_i ** 3) * y[j]
        res1[2, i] = v * x_i
        res1[3, i] = res1[1, i] - res1[2, i]

    print("Решение интегрального уравнения Фредгольма квадратурным методом")
    for row in res1:
        print(" ".join(f"{value:14.5f}" for value in row))

if __name__ == "__main__":
    main()