---

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid, simpson, quad

# Параметры задачи
a = 0.5
b = 1.5
n = 10
h = (b - a) / n  # 0.1

# Подынтегральная функция
def f(x):
    return 6 * x / np.sqrt(0.1 * x**2 + 1.8)

# 1. Левые прямоугольники (вперёд)
x_left = np.arange(a, b, h)          # от a до b-h
S_left = h * np.sum(f(x_left))

# 2. Правые прямоугольники (назад)
x_right = np.arange(a + h, b + h, h) # от a+h до b
S_right = h * np.sum(f(x_right))

# 3. Средние прямоугольники
x_mid = np.arange(a + h/2, b - h/2 + h/10, h)  # от a+h/2 до b-h/2
# +h/10 для компенсации погрешности округления (b-h/2 = 1.45)
S_mid = h * np.sum(f(x_mid))

# 4. Метод трапеций (ручная реализация)
x_trap = np.arange(a, b, h)          # от a до b-h
S_trap = h * np.sum((f(x_trap) + f(x_trap + h)) / 2)

# 5. Метод Симпсона (ручная реализация)
x_simp2 = np.arange(a + h, b, h)     # точки для коэффициента 2
x_simp4 = np.arange(a + h/2, b, h)   # точки для коэффициента 4 (до b-h/2)
S_simp = h/6 * (f(a) + f(b) + 2 * np.sum(f(x_simp2)) + 4 * np.sum(f(x_simp4)))

# 6. Встроенный метод трапеций (scipy / numpy)
x_all = np.linspace(a, b, n+1)       # равномерная сетка из 11 точек
y_all = f(x_all)
S_trapz = trapezoid(y_all, x_all)    # или np.trapz(y_all, x_all)

# 7. Встроенный метод Симпсона
S_simpson = simpson(y_all, x=x_all)

# 8. Адаптивный квадратурный метод quad
S_quad, _ = quad(f, a, b)

# Вывод результатов
print("Результаты вычисления интеграла ∫_{0.5}^{1.5} (6x / √(0.1x²+1.8)) dx")
print(f"1. Левые прямоугольники:   {S_left:.8f}")
print(f"2. Правые прямоугольники:  {S_right:.8f}")
print(f"3. Средние прямоугольники: {S_mid:.8f}")
print(f"4. Трапеции (ручной):      {S_trap:.8f}")
print(f"5. Симпсон (ручной):       {S_simp:.8f}")
print(f"6. scipy.trapezoid:        {S_trapz:.8f}")
print(f"7. scipy.simpson:          {S_simpson:.8f}")
print(f"8. scipy.quad:             {S_quad:.8f}")

# Построение графика функции
x_fine = np.linspace(a, b, 200)
y_fine = f(x_fine)

plt.figure(figsize=(10, 5))
plt.plot(x_fine, y_fine, 'b-', linewidth=2, label='f(x) = 6x / √(0.1x²+1.8)')
plt.fill_between(x_fine, y_fine, alpha=0.2, color='blue')
plt.title('График подынтегральной функции на отрезке [0.5, 1.5]')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.legend()
plt.tight_layout()
plt.show()