---

import numpy as np
import matplotlib.pyplot as plt


def solve_bvp(h, a=0.4, b=0.7, bc_system=2):

    n = int(round((b - a) / h))
    x = np.linspace(a, b, n + 1)

    p = -3.0                        
    q = lambda xk: 1.0 / xk         
    f = 1.0                         

    A = np.zeros((n + 1, n + 1))
    rhs = np.zeros(n + 1)

    A[0, 0] = 1.0
    rhs[0] = 2.0

    for k in range(1, n):
        xk = x[k]
        A[k, k-1] = 1.0 / h**2 - p / (2.0 * h)
        A[k, k]   = -2.0 / h**2 + q(xk)
        A[k, k+1] = 1.0 / h**2 + p / (2.0 * h)
        rhs[k] = f

    if bc_system == 1:

        A[n, n]   = h + 2.0
        A[n, n-1] = -2.0
        rhs[n] = 0.7 * h
    else:

        A[n, n]   = h + 3.0
        A[n, n-1] = -4.0
        A[n, n-2] = 1.0
        rhs[n] = 0.7 * h

    y = np.linalg.solve(A, rhs)
    return x, y



h_coarse = 0.05
h_fine   = 0.025

x1c, y1c = solve_bvp(h_coarse, bc_system=1)
x1f, y1f = solve_bvp(h_fine,   bc_system=1)

x2c, y2c = solve_bvp(h_coarse, bc_system=2)
x2f, y2f = solve_bvp(h_fine,   bc_system=2)

def max_diff(xc, yc, xf, yf):
    m = 0.0
    for xi, yi in zip(xc, yc):
        idx = np.argmin(np.abs(xf - xi))
        m = max(m, abs(yi - yf[idx]))
    return m

eps1 = max_diff(x1c, y1c, x1f, y1f)
eps2 = max_diff(x2c, y2c, x2f, y2f)


diff_coarse = max(abs(ya - yb) for ya, yb in zip(y1c, y2c))
diff_fine   = max(abs(ya - yb) for ya, yb in zip(y1f, y2f))

def print_solution(title, x, y):
   
    print(title)
    
    print(f"{'k':>3} {'x_k':>8} {'y_k':>16}")
   
    for k in range(len(x)):
        print(f"{k:>3} {x[k]:>8.4f} {y[k]:>16.10f}")
    print()


print_solution(f"Система 1 (O(h)),  h = {h_coarse}", x1c, y1c)
print_solution(f"Система 1 (O(h)), h = {h_fine}",   x2f, y2f)
print_solution(f"Система 2 (O(h²)), h = {h_coarse}", x2c, y2c)
print_solution(f"Система 2 (O(h²)), h = {h_fine}",   x2f, y2f)


print(" Система 1 при h = 0.05 и h = 0.025")

print(f"{'x':>8} {'y(h=0.05)':>16} {'y(h=0.025)':>16} {'|разность|':>14}")

for xi, ya, yb in zip(x1c, y1c, y1f):
    print(f"{xi:>8.4f} {ya:>16.10f} {yb:>16.10f} {abs(ya - yb):>14.3e}")
print(f"\nε₁ ≈ {eps1:.6e}\n")

print("Система 2 при h = 0.05 и h = 0.025")
print(f"{'x':>8} {'y(h=0.05)':>16} {'y(h=0.025)':>16} {'|разность|':>14}")
for xi, ya, yb in zip(x2c, y2c, x2f):
    print(f"{xi:>8.4f} {ya:>16.10f} {yb:>16.10f} {abs(ya - yb):>14.3e}")
print(f"\nε₂ ≈ {eps2:.6e}\n")


print("Сравнение Системы 1 и Системы 2 при h = 0.05")
print(f"{'x':>8} {'Система 1':>16} {'Система 2':>16} {'|разность|':>14}")

for xi, ya, yb in zip(x1c, y1c, y2c):   
    print(f"{xi:>8.4f} {ya:>16.10f} {yb:>16.10f} {abs(ya-yb):>14.3e}")
print()

def max_diff(xc, yc, xf, yf):
    m = 0.0
    for xi, yi in zip(xc, yc):
        idx = np.argmin(np.abs(xf - xi))
        m = max(m, abs(yi - yf[idx]))
    return m

eps1 = max_diff(x1c, y1c, x1f, y1f)
eps2 = max_diff(x2c, y2c, x2f, y2f)



print("Сравнение Системы 1 и Системы 2 при h = 0.025")
print(f"{'x':>8} {'Система 1':>16} {'Система 2':>16} {'|разность|':>14}")



for xi, ya, yb in zip(x1f, y1f, y2f):
    print(f"{xi:>8.4f} {ya:>16.10f} {yb:>16.10f} {abs(ya - yb):>14.3e}")
print()

diff_coarse = max(abs(ya - yb) for ya, yb in zip(y1c, y2c))
diff_fine   = max(abs(ya - yb) for ya, yb in zip(y1f, y2f))


print("Как меняется разница между Системой 1 и Системой 2")
print("при уменьшении шага")

print(f"max|С1 - С2| при h = 0.05 : {diff_coarse:.6e}")
print(f"max|С1 - С2| при h = 0.025: {diff_fine:.6e}")




print("Оценка точности (сравнение h = 0.05 и h = 0.025)")

print(f"Система 1 (O(h))  : ε ≈ {eps1:.6e}")
print(f"Система 2 (O(h²)) : ε ≈ {eps2:.6e}")
print()

if eps2 > 0:
    order = np.log2(eps1 / eps2)
    print(f"Наблюдаемый порядок сходимости Системы 2: ~{order:.2f}")
    print()


fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

axes[0].plot(x1c, y1c, 'o-',  label='Система 1 (O(h))',  markersize=7)
axes[0].plot(x2c, y2c, 's--', label='Система 2 (O(h²))', markersize=7)
axes[0].set_xlabel('x')
axes[0].set_ylabel('y(x)')
axes[0].set_title('Сравнение двух систем аппроксимации ГУ (h = 0.05)')
axes[0].grid(True)
axes[0].legend()

axes[1].plot(x2c, y2c, 'o-',  label=f'h = {h_coarse}', markersize=7)
axes[1].plot(x2f, y2f, 's--', label=f'h = {h_fine}',   markersize=4)
axes[1].set_xlabel('x')
axes[1].set_ylabel('y(x)')
axes[1].set_title('Сходимость Системы 2 при уменьшении шага')
axes[1].grid(True)
axes[1].legend()

fig2 = plt.figure(2, figsize=(9, 6))

ax3 = fig2.add_subplot(1, 1, 1)     
ax3.plot(x1c, y1c, 'o-',  color='tab:blue',
         label=f'Система 1, h = {h_coarse}', markersize=8)
ax3.plot(x1f, y1f, 's--', color='tab:red',
         label=f'Система 1, h = {h_fine}',   markersize=5)
ax3.set_xlabel('x')
ax3.set_ylabel('y(x)')
ax3.set_title('Система 1 (O(h)): решения при двух шагах')
ax3.grid(True)
ax3.legend()

fig2.tight_layout()
fig2.savefig('lab4_figure2.png', dpi=150)

plt.tight_layout()
plt.savefig('system1_plots.png', dpi=150)
plt.show()

plt.tight_layout()
plt.savefig('lab4_variant2.png', dpi=150)
plt.show()