---

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spi

# ------------------------------------------------------------
# Параметры задачи (можно легко изменить)
# ------------------------------------------------------------
def f(x):
    """Подынтегральная функция. Замените на свою при необходимости."""
    return x**2

a = 0.5          # начало отрезка
b = 1.5          # конец отрезка
n = 10           # число разбиений (чётное для Симпсона)
h = (b - a) / n  # шаг

# ------------------------------------------------------------
# 1. Левые прямоугольники
# ------------------------------------------------------------
def left_rect(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b - h, n)   # n точек от a до b-h
    return h * np.sum(f(x))

S_left = left_rect(f, a, b, n)

# ------------------------------------------------------------
# 2. Правые прямоугольники
# ------------------------------------------------------------
def right_rect(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a + h, b, n)   # n точек от a+h до b
    return h * np.sum(f(x))

S_right = right_rect(f, a, b, n)

# ------------------------------------------------------------
# 3. Средние прямоугольники
# ------------------------------------------------------------
def mid_rect(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a + h/2, b - h/2, n)
    return h * np.sum(f(x))

S_mid = mid_rect(f, a, b, n)

# ------------------------------------------------------------
# 4. Метод трапеций (собственная реализация)
# ------------------------------------------------------------
def trapezoid(f, a, b, n):
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    return h * (np.sum(y) - 0.5 * (y[0] + y[-1]))

S_trap = trapezoid(f, a, b, n)

# ------------------------------------------------------------
# 5. Метод Симпсона (собственная реализация)
# ------------------------------------------------------------
def simpson(f, a, b, n):
    if n % 2 != 0:
        raise ValueError("Число интервалов n должно быть чётным")
    h = (b - a) / n
    x = np.linspace(a, b, n + 1)
    y = f(x)
    # коэффициенты: 1, 4, 2, 4, 2, ..., 4, 1
    coeff = np.ones(n + 1)
    coeff[1:-1:2] = 4   # нечётные индексы внутри
    coeff[2:-1:2] = 2   # чётные индексы внутри
    return h / 3 * np.sum(coeff * y)

S_simp = simpson(f, a, b, n)

# ------------------------------------------------------------
# 6. Библиотечный метод трапеций (scipy)
# ------------------------------------------------------------
x_dense = np.linspace(a, b, 1000)   # для точности возьмём много точек
y_dense = f(x_dense)
S_trap_sp = spi.trapezoid(y_dense, x_dense)

# ------------------------------------------------------------
# 7. Библиотечный метод Симпсона (scipy)
# ------------------------------------------------------------
S_simp_sp = spi.simpson(y_dense, x=x_dense)

# ------------------------------------------------------------
# 8. Точное значение (quad)
# ------------------------------------------------------------
S_exact, err = spi.quad(f, a, b)

# ------------------------------------------------------------
# Вывод результатов
# ------------------------------------------------------------
print("Интеграл f(x) = x^2 на отрезке [0.5, 1.5]")
print("=" * 50)
print(f"Левые прямоугольники      : {S_left:.8f}")
print(f"Правые прямоугольники     : {S_right:.8f}")
print(f"Средние прямоугольники    : {S_mid:.8f}")
print(f"Трапеции (своя реализация): {S_trap:.8f}")
print(f"Симпсон (своя реализация) : {S_simp:.8f}")
print(f"Трапеции (scipy)          : {S_trap_sp:.8f}")
print(f"Симпсон (scipy)           : {S_simp_sp:.8f}")
print(f"Точное значение (quad)    : {S_exact:.8f}")
print("=" * 50)

# ------------------------------------------------------------
# Построение графика
# ------------------------------------------------------------
x_plot = np.linspace(a, b, 200)
y_plot = f(x_plot)

plt.figure(figsize=(8, 5))
plt.plot(x_plot, y_plot, 'b-', linewidth=2, label=r'$f(x) = x^2$')
plt.fill_between(x_plot, y_plot, alpha=0.2)
plt.title('График функции на интервале [0.5, 1.5]')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend()
plt.tight_layout()
plt.show()