№10 Метод неопределенных коэффициентов решения краевой задачи для ОДУ 2-го порядка

def f(x):
    V = 8
    return 4 * pow(x, 4) - 3 * V * pow(x, 3) + 6 * x - 2 * V

def p(x):
    return x * x

def q(x):
    return x

def check_y(x):
    V = 8
    return x * x * (x - V)

def phi(x, k):
    k += 1
    T = 8
    return pow(x, k + 1) * (x - T)

def phi1(x, k):
    k += 1
    T = 8
    return pow(x, k + 1) * (k + 2) - pow(x, k) * (k + 1) * T

def phi2(x, k):
    k += 1
    T = 8
    return pow(x, k) * (k + 2) * (k + 1) - pow(x, k - 1) * (k + 1) * k * T

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]) < 1e-6:
            print("Division by zero!")
            return []

        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][mx], sole[i][j] = sole[i][j], sole[i][mx]

        if abs(sole[j][mx]) < 1e-6:
            print("Division by zero!")
            return []

        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_rhs(A, B):
    for i in range(len(A)):
        A[i].append(B[i])
    return gauss(A)

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

        for j in range(i, m):
            print(f"{vy[j]:.5f}\t", end="")
        print()

        for j in range(i, m):
            print(f"{check_y(vx[j]):.5f}\t", end="")
        print()

        for j in range(i, m):
            print(f"{vy[j] - check_y(vx[j]):.5f}\t", end="")
        print()

        print()

def solve():
    n = 5
    l, r = 0, 8  # T = V = 16
    h = (r - l) / (n + 1)
    vx = [l + i * h for i in range(1, n + 1)]
    A = [[0] * n for _ in range(n)]
    B = [0] * n
    for j in range(n):
        x = vx[j]
        for k in range(n):
            A[j][k] = phi2(x, k) + p(x) * phi1(x, k) + q(x) * phi(x, k)
        B[j] = f(x)

    va = gauss_with_rhs(A, B)
    vx.insert(0, 0)
    vx.append(8)
    vy = [0] * len(vx)
    for i in range(len(vx)):
        for k in range(n):
            vy[i] += va[k] * phi(vx[i], k)

    print_result(vx, vy)

if __name__ == "__main__":
    solve()