https://pastein.ru/t/XrX

  скопируйте уникальную ссылку для отправки

Загрузка данных


from math import sin, cos, sqrt, pi, inf, log
import numpy as np

# A, B = 6, 9
A, B = 0, 1
EPS = 10**-7


def f1(x):
    # return sqrt((9-2*x) / (2*x-21))
    return log(1+x) / ((3*x+2)**2)


# метод трапеций
def trapezoidal_rule(f, a, b, n):
    h = (b - a) / n
    s = 0.5 * (f(a) + f(b))
    for i in range(1, n):
        s += f(a + i * h)
    return s * h


# метод Симпсона
def simpson_rule(f, a, b, n):
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n, 2):
        s += 4 * f(a + i * h)
    for i in range(2, n-1, 2):
        s += 2 * f(a + i * h)
    return s * h / 3


def runge_eps(s1, s2, m):
    return abs(s2 - s1) / (2**m - 1)


def print_row(n, h, i, r):
    add_tab = '\t' if n < 10 else ''
    print(f'N={n}{add_tab}\th={h:.12f}\tI={i:.12f}\tR={r:.12f}')


# метод трапеций с точностью eps
def trapezoidal_rule_eps(f, a, b, eps):
    print('trapezoidal_rule_eps')
    n = 2
    m = 2  # порядок точности
    s1 = inf
    s2 = trapezoidal_rule(f, a, b, 1)
    while runge_eps(s1, s2, m) > eps:
        s1 = s2
        s2 = trapezoidal_rule(f, a, b, n)
        print_row(n, (b-a)/n, s2, runge_eps(s1, s2, m))
        n *= 2
    print()
    return s2, n // 2, runge_eps(s1, s2, m)


# метод Симпсона с точностью eps
def simpson_rule_eps(f, a, b, eps):
    print('simpson_rule_eps')
    n = 2
    m = 4  # порядок точности
    s1, s2 = inf, simpson_rule(f, a, b, 1)
    while runge_eps(s1, s2, m) > eps:
        s1, s2 = s2, simpson_rule(f, a, b, n)
        print_row(n, (b-a)/n, s2, runge_eps(s1, s2, m))
        n *= 2
    print()
    return s2, n // 2, runge_eps(s1, s2, m)


# K = 4
# gauss_coeff = [
#     [sqrt(3/7-2/7*sqrt(6/5)), (18+sqrt(30))/36],
#     [-sqrt(3/7-2/7*sqrt(6/5)), (18+sqrt(30))/36],
#     [sqrt(3/7+2/7*sqrt(6/5)), (18-sqrt(30))/36],
#     [-sqrt(3/7+2/7*sqrt(6/5)), (18-sqrt(30))/36],
# ]
K = 5
gauss_coeff = [
    [0, 128/255],
    [1/3 * sqrt(5-2*sqrt(10/7)), (322+13*sqrt(70))/900],
    [-1/3 * sqrt(5-2*sqrt(10/7)), (322+13*sqrt(70))/900],
    [1/3 * sqrt(5+2*sqrt(10/7)), (322-13*sqrt(70))/900],
    [-1/3 * sqrt(5+2*sqrt(10/7)), (322-13*sqrt(70))/900]
]
def nast_rule(f, a, b):
    r = 0
    for x_g, w_g in gauss_coeff:
        # конвертируем x на наш отрезок
        x = (a + b) / 2.0 + (b - a) / 2.0 * x_g
        r += f(x) * w_g
    r *= (b - a) / 2
    return r


def main():
    # prec_value = pi
    prec_value = 1/15*(5*log(5)-11*log(2))
    trapez_res = trapezoidal_rule(f1, A, B, 1000)
    simpson_res = simpson_rule(f1, A, B, 1000)
    trapez_eps_res, trapez_steps, trapez_eps_runge = trapezoidal_rule_eps(f1, A, B, EPS)
    simpson_eps_res, simpson_steps, simpson_eps_runge = simpson_rule_eps(f1, A, B, EPS)

    print(f"Точное значение: {prec_value}")
    print()

    print(f"Метод трапеций с 1000 разбиениями: {trapez_res}")
    print(f"Реальная погрешность: {abs(prec_value - trapez_res):.12f}")
    print(f"Метод Симпсона с 1000 разбиениями: {simpson_res}")
    print(f"Реальная погрешность: {abs(prec_value - simpson_res):.12f}")
    print()

    print(f"Метод трапеций с точностью {EPS}: {trapez_eps_res}")
    print(f"Количество разбиений: {trapez_steps}")
    print(f"Погрешность по Рунге: {trapez_eps_runge:.12f}")
    print(f"Реальная погрешность: {abs(prec_value - trapez_eps_res):.12f}")
    print(f"Метод Симпсона с точностью {EPS}: {simpson_eps_res}")
    print(f"Количество разбиений: {simpson_steps}")
    print(f"Погрешность по Рунге: {simpson_eps_runge:.12f}")
    print(f"Реальная погрешность: {abs(prec_value - simpson_eps_res):.12f}")
    print()

    nast_res = nast_rule(f1, A, B)
    print(f"КФ НАСТ с k={K}: {nast_res}")
    print(f"Реальная погрешность: {abs(prec_value - nast_res):.12f}")


if __name__ == '__main__':
    main()