import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import time

eps = 0.9468
L   = 30001
N   = 14
x0  = np.array([1.0, 0.0])

t  = np.linspace(-eps, eps, L)
h  = t[1] - t[0]
i0 = L // 2
inicio = time.perf_counter()

def F(t, x):
    dx1 = x[:, 1]
    dx2 = -x[:, 0] + 2.0 * np.sin(t)
    return np.column_stack([dx1, dx2])

def x_exata(t):
    x1 = (1.0 - t) * np.cos(t) + np.sin(t)
    x2 = -(1.0 - t) * np.sin(t)
    return np.column_stack([x1, x2])

# Iterados de Picard

y_tp = np.zeros((N+1, L, 2))
y_tp[0, :, :] = x0

y_rt = np.zeros((N+1, L, 2))
y_rt[0, :, :] = x0

for n in range(N):
    f_tp = F(t, y_tp[n])
    f_rt = F(t, y_rt[n])

    I_tp_pos = np.zeros((L, 2))
    I_rt_pos = np.zeros((L, 2))
    for k in range(i0 + 1, L):
        I_tp_pos[k] = I_tp_pos[k-1] + (f_tp[k-1] + f_tp[k]) / 2.0 * h
        I_rt_pos[k] = I_rt_pos[k-1] +  f_rt[k-1]             * h

    I_tp_neg = np.zeros((L, 2))
    I_rt_neg = np.zeros((L, 2))
    for k in range(i0 - 1, -1, -1):
        I_tp_neg[k] = I_tp_neg[k+1] + (f_tp[k] + f_tp[k+1]) / 2.0 * h
        I_rt_neg[k] = I_rt_neg[k+1] +  f_rt[k]               * h

    y_tp[n+1, i0, :]  = x0
    y_tp[n+1, i0+1:]  = x0 + I_tp_pos[i0+1:]
    y_tp[n+1, :i0]    = x0 - I_tp_neg[:i0]

    y_rt[n+1, i0, :]  = x0
    y_rt[n+1, i0+1:]  = x0 + I_rt_pos[i0+1:]
    y_rt[n+1, :i0]    = x0 - I_rt_neg[:i0]


# Tabela de erros máximos

x_ref    = x_exata(t)
erros_tp = np.zeros(N+1)
erros_rt = np.zeros(N+1)

print(f"{'n':<4} {'trapézio':>12} {'retângulo':>12}")
print("-" * 30)
for n in range(N+1):
    erros_tp[n] = np.max(np.abs(y_tp[n] - x_ref))
    erros_rt[n] = np.max(np.abs(y_rt[n] - x_ref))
    print(f"{n:<4} {erros_tp[n]:>12.2e} {erros_rt[n]:>12.2e}")

cores = cm.viridis(np.linspace(0, 1, N+1))
sm = cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=2, vmax=N))

#  Figura 1: Trapézio (x1 e x2 sobrepostos num único eixo) 
fig1, ax1 = plt.subplots(figsize=(8, 5))
for n in range(2, N+1):
    ax1.plot(t, y_tp[n, :, 0], color=cores[n], lw=0.8)
    ax1.plot(t, y_tp[n, :, 1], color=cores[n], lw=0.8, linestyle=':')

ax1.plot(t, x_ref[:, 0], linestyle='--', color='red',    linewidth=2, label='exata $x_1$')
ax1.plot(t, x_ref[:, 1], linestyle='--', color='darkred', linewidth=2, label='exata $x_2$')

ax1.plot([], [], color='gray', lw=1.2,                  label='iterados $x_1$ (sólido)')
ax1.plot([], [], color='gray', lw=1.2, linestyle=':',   label='iterados $x_2$ (pontilhado)')
ax1.set_title("Trapézio")
ax1.set_xlabel("$t$"); ax1.set_ylabel("$x$")
ax1.legend(fontsize=8); ax1.grid(True, alpha=0.4)
plt.colorbar(cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=2, vmax=N)),
             ax=ax1, label='iterado $n$')
plt.tight_layout()

plt.show()

#  Figura 2: Retângulo (x1 e x2 sobrepostos num único eixo) 
fig2, ax2 = plt.subplots(figsize=(8, 5))
for n in range(2, N+1):
    ax2.plot(t, y_rt[n, :, 0], color=cores[n], lw=0.8)
    ax2.plot(t, y_rt[n, :, 1], color=cores[n], lw=0.8, linestyle=':')
ax2.plot(t, x_ref[:, 0], linestyle='--', color='red',    linewidth=2, label='exata $x_1$')
ax2.plot(t, x_ref[:, 1], linestyle='--', color='darkred', linewidth=2, label='exata $x_2$')
ax2.plot([], [], color='gray', lw=1.2,                  label='iterados $x_1$ (sólido)')
ax2.plot([], [], color='gray', lw=1.2, linestyle=':',   label='iterados $x_2$ (pontilhado)')
ax2.set_title("Retângulo")
ax2.set_xlabel("$t$"); ax2.set_ylabel("$x$")
ax2.legend(fontsize=8); ax2.grid(True, alpha=0.4)
plt.colorbar(cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=2, vmax=N)),
             ax=ax2, label='iterado $n$')
plt.tight_layout()
plt.show()

#  Figura 3: Erro pontual — Trapézio 
fig3, ax3 = plt.subplots(figsize=(7, 5))
ax3.plot(t, np.abs(y_tp[N, :, 0] - x_ref[:, 0]), label="trapézio $x_1$", color='steelblue')
ax3.plot(t, np.abs(y_tp[N, :, 1] - x_ref[:, 1]), label="trapézio $x_2$", color='tomato')
ax3.set_yscale('log')
ax3.set_xlabel("$t$"); ax3.set_ylabel("erro pontual")
ax3.set_title(f"Erro pontual — Trapézio, iterado $n = {N}$")
ax3.legend(); ax3.grid(True, alpha=0.4)
plt.tight_layout()
plt.show()

#  Figura 4: Erro pontual — Retângulo 
fig4, ax4 = plt.subplots(figsize=(7, 5))
ax4.plot(t, np.abs(y_rt[N, :, 0] - x_ref[:, 0]), label="retângulo $x_1$", color='steelblue')
ax4.plot(t, np.abs(y_rt[N, :, 1] - x_ref[:, 1]), label="retângulo $x_2$", color='tomato')
ax4.set_yscale('log')
ax4.set_xlabel("$t$"); ax4.set_ylabel("erro pontual")
ax4.set_title(f"Erro pontual — Retângulo, iterado $n = {N}$")
ax4.legend(); ax4.grid(True, alpha=0.4)
plt.tight_layout()
plt.show()


#  Tabela 9 (EDO5): tolerâncias avaliadas com δ* ótimo 
TOLS_TABELA9 = [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12]

#  Figura Raio de convergência: intervalo de δ para as curvas 
DELTA_MIN_CURVA = 0.01
DELTA_MAX_CURVA = 5.0

#  Figura Iterados necessários: tolerâncias das curvas e limite do eixo x 
TOLS_CURVAS    = [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12]
DELTA_MAX_PLOT = 5.0

#  Figura Barras / Tabela comparativa: tolerâncias 
TOLS_BARRAS = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]

delta_curva = np.linspace(DELTA_MIN_CURVA, DELTA_MAX_CURVA, 3000)
C_vals      = np.ones_like(delta_curva)
M_vals      = np.sqrt(delta_curva**2 +
                      (delta_curva - 1.0 + 2.0 * np.sin(delta_curva))**2)
eps_ratio   = delta_curva / M_vals
eps_vals    = np.minimum(delta_curva, eps_ratio)

idx_opt = np.argmax(eps_vals)
d_opt   = delta_curva[idx_opt]
eps_opt = eps_vals[idx_opt]
M_opt   = M_vals[idx_opt]
C_opt   = C_vals[idx_opt]

print("\n--- Análise teórica de convergência (sistema x'' = -x + 2 sen t) ---")
print(f"  C(δ) constante = {C_opt:.6f}  (exato: 1.0)")
print(f"  δ* ótimo       = {d_opt:.4f}")
print(f"  ε(δ*)          = {eps_opt:.4f}")
print(f"  M(δ*)          = {M_opt:.4f}")

def nmax_scalar(C_val, eps_val, M_val, tol):
    """Menor n tal que M·(C·ε)^n/n! < tol (converge sempre pelo crescimento de n!)."""
    termo = float(M_val)
    n = 0
    while termo >= tol and n < 10000:
        n += 1
        termo *= (C_val * eps_val) / n
    return n

# ── Tabela 9 — N_max(τ) com δ* ótimo ────────────────────────────────────────
print(f"\n--- Tabela 9 (EDO5): N_max(τ) com δ* = {d_opt:.4f}, ε* = {eps_opt:.4f} ---")
print(f"{'Tolerância τ':>15} {'N_max':>8}")
print("-" * 26)
for tol in TOLS_TABELA9:
    nm = nmax_scalar(C_opt, eps_opt, M_opt, tol)
    print(f"{tol:>15.0e} {nm:>8}")

# ── Tabela 13 — N_max comparativo EDO4 vs EDO5 ──────────────────────────────
# EDO4: δ*=1, ε*=0.25, M*=4, C*=4  (F=x², analítico — já calculado no EDO4)
C4, e4, M4 = 4.0, 0.25, 4.0
print(f"\n--- Tabela 13: N_max comparativo EDO4 vs EDO5 ---")
print(f"{'Tolerância τ':>15} {'N_max EDO4':>12} {'N_max EDO5 (2D)':>16}")
print("-" * 46)
for tol in TOLS_TABELA9:
    nm4 = nmax_scalar(C4, e4, M4, tol)
    nm5 = nmax_scalar(C_opt, eps_opt, M_opt, tol)
    print(f"{tol:>15.0e} {nm4:>12} {nm5:>16}")

#  Figura 5: Raio de convergência 
fig5, ax5 = plt.subplots(figsize=(7, 5))
ax5.plot(delta_curva, eps_vals, label=r"$\varepsilon(\delta)$",  color='steelblue', lw=2)
ax5.plot(delta_curva, M_vals,   label=r"$M(\delta)$",            color='tomato',    lw=2)
ax5.plot(delta_curva, C_vals,   label=r"$C(\delta)=1$",          color='seagreen',  lw=2, linestyle='--')
ax5.axvline(d_opt,   color='k',         linestyle='--', lw=1.2, label=rf"$\delta^*={d_opt:.3f}$")
ax5.axhline(eps_opt, color='steelblue', linestyle=':',  lw=1.2, label=rf"$\varepsilon^*={eps_opt:.3f}$")
ax5.set_xlabel(r"$\delta$"); ax5.set_ylabel("")
ax5.set_xlim(0, DELTA_MAX_CURVA); ax5.set_ylim(0, 5)
ax5.set_title(r"Raio de convergência $\varepsilon(\delta)$, $M(\delta)$, $C(\delta)$")
ax5.legend(fontsize=9); ax5.grid(True, alpha=0.4)
plt.tight_layout()
plt.show()

#  Figura 6: Iterados necessários N_max(δ, τ) 
cores_tol = plt.cm.plasma(np.linspace(0.1, 0.9, len(TOLS_CURVAS)))
mask_plot = delta_curva <= DELTA_MAX_PLOT
d_plot    = delta_curva[mask_plot]

fig6, ax6 = plt.subplots(figsize=(7, 5))
for tol, cor in zip(TOLS_CURVAS, cores_tol):
    Nmax_arr = np.array([nmax_scalar(C_vals[i], eps_vals[i], M_vals[i], tol)
                         for i in range(len(d_plot))])
    ax6.plot(d_plot, Nmax_arr, color=cor, lw=1.8, label=f"$\\tau={tol:.0e}$")
ax6.axvline(d_opt, color='k', linestyle='--', lw=1.2, label=rf"$\delta^*={d_opt:.3f}$")
ax6.set_xlabel(r"$\delta$"); ax6.set_ylabel(r"$N_{\max}(\tau)$")
ax6.set_ylim(0, 30); ax6.set_xlim(0, DELTA_MAX_PLOT)
ax6.set_title(r"Iterados necessários $N_{\max}(\delta,\tau)$")
ax6.legend(fontsize=8); ax6.grid(True, alpha=0.4)
plt.tight_layout()
plt.show()

#  Figura 7: Convergência (erro máx vs n) + N_max teórico vs numérico 
n_vals   = np.arange(N+1)
cota_teo = np.zeros(N+1)
cota_teo[0] = M_opt
for n in range(1, N+1):
    cota_teo[n] = cota_teo[n-1] * (C_opt * eps_opt) / n

tols_cmp = TOLS_BARRAS
n_teo    = []
n_num_tp = []
n_num_rt = []

for tol in tols_cmp:
    n_teo.append(nmax_scalar(C_opt, eps_opt, M_opt, tol))
    cruzou_tp = next((n for n in range(N+1)
                      if np.max(np.abs(y_tp[n] - x_ref)) < tol), None)
    cruzou_rt = next((n for n in range(N+1)
                      if np.max(np.abs(y_rt[n] - x_ref)) < tol), None)
    n_num_tp.append(cruzou_tp if cruzou_tp is not None else N+1)
    n_num_rt.append(cruzou_rt if cruzou_rt is not None else N+1)

print(f"\n{'τ':>10} {'N_teo':>8} {'n_tp':>10} {'n_rt':>10}")
print("-" * 42)
for tol, nt, ntp, nrt in zip(tols_cmp, n_teo, n_num_tp, n_num_rt):
    flag_tp = "†" if ntp > N else ""
    flag_rt = "†" if nrt > N else ""
    print(f"{tol:>10.0e} {nt:>8} {str(ntp)+flag_tp:>10} {str(nrt)+flag_rt:>10}")
print("  † não atingiu a tolerância em N iterados")

#    Figura 7: Convergência dos iterados (erro máx vs n) 
fig7, ax7 = plt.subplots(figsize=(8, 5))
ax7.semilogy(n_vals, erros_tp, 'o-',  color='steelblue', label="trapézio (numérico)",  lw=1.8)
ax7.semilogy(n_vals, erros_rt, 's--', color='tomato',    label="retângulo (numérico)", lw=1.8)
ax7.semilogy(n_vals, cota_teo, '^:',  color='seagreen',
             label=rf"cota teórica ($\delta^*={d_opt:.2f}$, $\varepsilon^*={eps_opt:.2f}$)", lw=1.8)
ax7.set_xlabel("iterado $n$"); ax7.set_ylabel("erro máximo")
ax7.set_title(r"Convergência dos iterados de Picard — $\varepsilon_{\rm num}=" + f"{eps}$")
ax7.legend(); ax7.grid(True, alpha=0.4, which='both')
plt.tight_layout()
plt.show()

#    Figura 8: N_max teórico vs cruzamento numérico (barras)    
fig8, ax8 = plt.subplots(figsize=(8, 5))
x_pos  = np.arange(len(tols_cmp))
w      = 0.25
labels = [f"$10^{{{int(np.log10(tol))}}}$" for tol in tols_cmp]
ax8.bar(x_pos - w, n_teo,    width=w, label=r"$N_{\max}$ teórico ($\delta^*$)", color='steelblue')
ax8.bar(x_pos,     n_num_tp, width=w, label="trapézio numérico",               color='tomato')
ax8.bar(x_pos + w, n_num_rt, width=w, label="retângulo numérico",              color='seagreen')
ax8.set_xticks(x_pos); ax8.set_xticklabels(labels)
ax8.set_xlabel(r"tolerância $\tau$"); ax8.set_ylabel("iterados necessários")
ax8.set_title(r"$x'' = -x + 2\sin t$,  $C(\delta)\equiv 1$")
ax8.legend(); ax8.grid(True, axis='y', alpha=0.4)
plt.tight_layout()
plt.show()

fim = time.perf_counter()
print(f"Tempo de execução: {fim - inicio:.2f} s")