---

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

# ========== Вариант 23 ==========
print("="*70)
print("Вариант 23: (x^2+1)y' + xy = x(x^2+1),  y(√2)=1,  x∈[√2, √2+1]")
print("="*70)

# Правая часть, приведённая к виду y' = f(x,y)
f23 = lambda x, y: x - (x/(x**2+1))*y

# Точное решение
y_true23 = lambda x: (x**2 + 1)/3

# Параметры
a23, b23 = np.sqrt(2), np.sqrt(2) + 1
y0_23 = 1.0
n_intervals = 5               # 5 разбиений -> 6 точек
x_grid23 = np.linspace(a23, b23, n_intervals+1)
h23 = (b23 - a23) / n_intervals

# --- Метод Эйлера ---
y_euler23 = np.zeros(n_intervals+1)
y_euler23[0] = y0_23
for i in range(n_intervals):
    y_euler23[i+1] = y_euler23[i] + h23 * f23(x_grid23[i], y_euler23[i])

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

# --- Решатель scipy (RK45) ---
sol23 = solve_ivp(f23, [a23, b23], [y0_23], method='RK45', t_eval=x_grid23)
y_solver23 = sol23.y[0]

# --- Точное решение ---
y_exact23 = y_true23(x_grid23)

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

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

# ========== Вариант 29 ==========
print("\n" + "="*70)
print("Вариант 29: x y y' = 1 - x^2,  y(0.5)=1.3,  x∈[0.5, 1.5]")
print("="*70)

# Правая часть: y' = (1 - x^2)/(x*y)
f29 = lambda x, y: (1 - x**2)/(x * y) if y != 0 else 0

# Точное решение (находим константу из начального условия)
x0_29, y0_29 = 0.5, 1.3
C29 = y0_29**2 - 2*np.log(x0_29) + x0_29**2   # из y² = 2 ln x - x² + C
y_true29 = lambda x: np.sqrt(2*np.log(x) - x**2 + C29)   # берём положительный корень

# Параметры
a29, b29 = 0.5, 1.5
n_intervals = 5
x_grid29 = np.linspace(a29, b29, n_intervals+1)
h29 = (b29 - a29) / n_intervals

# Метод Эйлера
y_euler29 = np.zeros(n_intervals+1)
y_euler29[0] = y0_29
for i in range(n_intervals):
    y_euler29[i+1] = y_euler29[i] + h29 * f29(x_grid29[i], y_euler29[i])

# Модифицированный метод Эйлера
y_mod29 = np.zeros(n_intervals+1)
y_mod29[0] = y0_29
for i in range(n_intervals):
    k1 = f29(x_grid29[i], y_mod29[i])
    y_pred = y_mod29[i] + h29 * k1
    k2 = f29(x_grid29[i] + h29, y_pred)
    y_mod29[i+1] = y_mod29[i] + h29 * (k1 + k2) / 2.0

# Решатель scipy
sol29 = solve_ivp(f29, [a29, b29], [y0_29], method='RK45', t_eval=x_grid29)
y_solver29 = sol29.y[0]

# Точное решение
y_exact29 = y_true29(x_grid29)

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

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