---

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

# ============================================================
# Постановка задачи (вариант 2)
#   y'' = x + 1,   y(4) = y(5) = 0
# ============================================================
a, b = 4.0, 5.0
L = b - a

f_rhs = lambda x: x + 1.0   # правая часть

# ============================================================
# Базисные функции: синусы (как в методичке 10.pdf)
#   phi_i(x) = sin(pi * i * (x-a) / (b-a))
# ============================================================
def phi(i, x):
    return np.sin(np.pi * i * (x - a) / L)

def phi_dd(i, x):
    # Вторая производная по x
    k = np.pi * i / L
    return -(k**2) * np.sin(k * (x - a))

N = 3   # число базисных функций

# ============================================================
# 1) МЕТОД КОЛЛОКАЦИЙ
# ============================================================
x_col = np.linspace(a, b, N+2)[1:-1]
A_col = np.zeros((N, N))
b_col = np.zeros(N)
for j, xj in enumerate(x_col):
    for i in range(1, N+1):
        A_col[j, i-1] = phi_dd(i, xj)
    b_col[j] = f_rhs(xj)
c_col = np.linalg.solve(A_col, b_col)

# ============================================================
# 2) МЕТОД НАИМЕНЬШИХ КВАДРАТОВ
# ============================================================
A_ls = np.zeros((N, N))
b_ls = np.zeros(N)
for i in range(1, N+1):
    for j in range(1, N+1):
        A_ls[i-1, j-1] = quad(lambda x: phi_dd(i, x) * phi_dd(j, x), a, b)[0]
    b_ls[i-1] = quad(lambda x: f_rhs(x) * phi_dd(i, x), a, b)[0]
c_ls = np.linalg.solve(A_ls, b_ls)

# ============================================================
# 3) МЕТОД ГАЛЁРКИНА
# ============================================================
A_gal = np.zeros((N, N))
b_gal = np.zeros(N)
for i in range(1, N+1):
    for j in range(1, N+1):
        A_gal[i-1, j-1] = quad(lambda x: phi_dd(i, x) * phi(j, x), a, b)[0]
    b_gal[i-1] = quad(lambda x: f_rhs(x) * phi(i, x), a, b)[0]
c_gal = np.linalg.solve(A_gal, b_gal)

# ============================================================
# Точное решение (для [4,5])
# ============================================================
y_exact = lambda x: x**3/6.0 + x**2/2.0 - (44.0/3.0)*x + 40.0

def y_approx(x, c):
    return sum(c[i-1] * phi(i, x) for i in range(1, N+1))

# ============================================================
# Вывод
# ============================================================
print("="*60)
print("  Коэффициенты (синусоидальный базис, N = 3)")
print("="*60)
print(f"  Коллокации : c = {np.round(c_col, 4)}")
print(f"  МНК        : c = {np.round(c_ls, 4)}")
print(f"  Галёркин   : c = {np.round(c_gal, 4)}")

xs = np.linspace(a, b, 11)
print("\n" + "="*85)
print(f"{'x':>6} | {'y_точн':>10} | {'Коллокации':>12} | {'МНК':>12} | {'Галёркин':>12} | {'max|err|':>10}")
print("-"*85)
err_max = 0.0
for x in xs:
    ye = y_exact(x)
    yc = y_approx(x, c_col)
    yl = y_approx(x, c_ls)
    yg = y_approx(x, c_gal)
    emax = max(abs(yc-ye), abs(yl-ye), abs(yg-ye))
    err_max = max(err_max, emax)
    print(f"{x:6.2f} | {ye:10.6f} | {yc:12.6f} | {yl:12.6f} | {yg:12.6f} | {emax:10.2e}")
print("="*85)
print(f"Максимальная погрешность: {err_max:.4e}")

# ============================================================
# График
# ============================================================
xx = np.linspace(a, b, 400)
plt.figure(figsize=(10, 6))
plt.plot(xx, y_exact(xx), 'k-',  lw=3, label='Точное решение')
plt.plot(xx, y_approx(xx, c_col), 'r--', lw=2, label='Коллокации')
plt.plot(xx, y_approx(xx, c_ls),  'b-.', lw=2, label='МНК')
plt.plot(xx, y_approx(xx, c_gal), 'g:',  lw=2.5, label='Галёркин')
plt.scatter(x_col, y_exact(x_col), c='r', s=80, zorder=5, label='Точки коллокации')

plt.xlabel('x'); plt.ylabel('y(x)')
plt.title('Вариант 2: y\'\' = x+1, y(4)=y(5)=0 (синусоидальный базис)')
plt.legend(); plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()