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()