---

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ========== Вариант 16 ==========
print("="*70)
print("Вариант 16: x y' - y = x^2 sin x,  y(0)=0,  x∈[π/2, π/2+1]")
print("="*70)

# Правая часть, приведённая к виду y' = f(x,y)
# Из xy' - y = x^2 sin x  =>  y' = y/x + x sin x
f16 = lambda x, y: y/x + x*np.sin(x) if x != 0 else 0

# Точное решение: y_top = -x cos x
y_true16 = lambda x: -x * np.cos(x)

# Параметры
a16, b16 = np.pi/2, np.pi/2 + 1
y0_16 = 0.0
n_intervals = 5                     # 5 разбиений -> 6 точек
x_grid16 = np.linspace(a16, b16, n_intervals+1)
h16 = (b16 - a16) / n_intervals

# --- Метод Эйлера ---
y_euler16 = np.zeros(n_intervals+1)
y_euler16[0] = y0_16
for i in range(n_intervals):
    y_euler16[i+1] = y_euler16[i] + h16 * f16(x_grid16[i], y_euler16[i])

# --- Модифицированный метод Эйлера (предиктор-корректор) ---
y_mod16 = np.zeros(n_intervals+1)
y_mod16[0] = y0_16
for i in range(n_intervals):
    k1 = f16(x_grid16[i], y_mod16[i])
    y_pred = y_mod16[i] + h16 * k1
    k2 = f16(x_grid16[i] + h16, y_pred)
    y_mod16[i+1] = y_mod16[i] + h16 * (k1 + k2) / 2.0

# --- Решатель scipy (RK45) ---
sol16 = solve_ivp(f16, [a16, b16], [y0_16], method='RK45', t_eval=x_grid16)
y_solver16 = sol16.y[0]

# --- Точное решение ---
y_exact16 = y_true16(x_grid16)

# Печать таблицы
print(f"{'x':>10} | {'Эйлер':>12} | {'Мод.Эйлер':>12} | {'solve_ivp':>12} | {'Точное':>12}")
print("-"*70)
for i in range(n_intervals+1):
    print(f"{x_grid16[i]:10.6f} | {y_euler16[i]:12.6f} | {y_mod16[i]:12.6f} | {y_solver16[i]:12.6f} | {y_exact16[i]:12.6f}")

# График
plt.figure(figsize=(10,5))
plt.plot(x_grid16, y_euler16, 'o-', label='Эйлер', markersize=4)
plt.plot(x_grid16, y_mod16, 's-', label='Модифицированный Эйлер', markersize=4)
plt.plot(x_grid16, y_solver16, '^-', label='solve_ivp (RK45)', markersize=4)
plt.plot(x_grid16, y_exact16, '*-', label='Точное решение', linewidth=2, markersize=5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Вариант 16: сравнение методов')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('variant16.png', dpi=150)
plt.show()