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

eps = 0.99
L   = 20001
N   = 30
x0  = 1.0

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

def F(t, x):
    return x**2

def x_exata(t):
    return 1.0 / (1.0 - t)


# Iterados de Picard — versão incremental 

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

y_rt = np.zeros((N+1, L))
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)
    I_rt_pos = np.zeros(L)
    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)
    I_rt_neg = np.zeros(L)
    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)

print(f"{'n':<4} {'trapézio':>12} {'retângulo':>12}")
print("-" * 30)
erros_tp = np.zeros(N+1)
erros_rt = np.zeros(N+1)
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))

#  Figura 1: Trapézio 
fig1, ax1 = plt.subplots(figsize=(7, 5))
for n in range(1, N+1):
    ax1.plot(t, y_tp[n], color=cores[n])
ax1.plot(t, x_ref, linestyle='--', color='red', linewidth=2, label='exata')
ax1.set_title("Trapézio")
ax1.set_xlabel("$t$"); ax1.set_ylabel("$x$")
ax1.set_xlim(0.7, 0.99); ax1.set_ylim(0, 30)
ax1.legend(); ax1.grid(True)
plt.colorbar(cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=1, vmax=N)),
             ax=ax1, label='iterado $n$')
plt.tight_layout()
plt.show()

#  Figura 2: Retângulo 
fig2, ax2 = plt.subplots(figsize=(7, 5))
for n in range(1, N+1):
    ax2.plot(t, y_rt[n], color=cores[n])
ax2.plot(t, x_ref, linestyle='--', color='red', linewidth=2, label='exata')
ax2.set_title("Retângulo")
ax2.set_xlabel("$t$"); ax2.set_ylabel("$x$")
ax2.set_xlim(0.7, 0.99); ax2.set_ylim(0, 30)
ax2.legend(); ax2.grid(True)
plt.colorbar(cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=1, vmax=N)),
             ax=ax2, label='iterado $n$')
plt.tight_layout()

plt.show()

#  Figura 3: Erro pontual 
fig3, ax3 = plt.subplots(figsize=(7, 5))
ax3.plot(t, np.abs(y_tp[N] - x_ref), label="erro trapézio")
ax3.plot(t, np.abs(y_rt[N] - x_ref), label="erro retângulo")
ax3.set_yscale('log')
ax3.set_xlabel("$t$"); ax3.set_ylabel("erro pontual")
ax3.set_title(f"Erro pontual — iterado $n = {N}$")
ax3.legend(); ax3.grid(True)
plt.tight_layout()

plt.show()


#  Tabela: δ avaliados 
DELTA_TABELA = [0.25, 0.50, 1.00, 1.50, 2.00]

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

#  intervalo de δ para as curvas 
DELTA_MIN_CURVA = 0.01
DELTA_MAX_CURVA = 3.0

#  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 = 2.0

# tolerâncias comparadas 
TOLS_BARRAS = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]



def M_func(d):   return (1.0 + d)**2
def C_func(d):   return 2.0 * (1.0 + d)
def eps_func(d): return d / M_func(d)

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

# δ* ótimo (máximo de ε(δ) = δ/(1+δ)²  →  δ* = 1 analiticamente)
delta_curva = np.linspace(DELTA_MIN_CURVA, DELTA_MAX_CURVA, 2000)
M_curva     = M_func(delta_curva)
C_curva     = C_func(delta_curva)
eps_curva   = eps_func(delta_curva)

idx_opt = np.argmax(eps_curva)
d_opt   = delta_curva[idx_opt]
eps_opt = eps_curva[idx_opt]
M_opt   = M_curva[idx_opt]
C_opt   = C_curva[idx_opt]

# ── Tabela 8 ─────────────────────────────────────────────────────────────────
print("\n--- Tabela 8: M(δ), C(δ), ε(δ), N_max(τ=1e-6) para F = x² ---")
print(f"{'δ':>6} {'M(δ)':>8} {'C(δ)':>8} {'ε(δ)':>8} {'N_max(τ=1e-6)':>15}")
print("-" * 52)
for d in DELTA_TABELA:
    M_d   = M_func(d)
    C_d   = C_func(d)
    e_d   = eps_func(d)
    nm    = nmax_scalar(C_d, e_d, M_d, 1e-6)
    print(f"{d:>6.2f} {M_d:>8.4f} {C_d:>8.4f} {e_d:>8.4f} {nm:>15}")

print(f"\n  δ* ótimo = {d_opt:.4f}  (analítico: 1.0000)")
print(f"  ε(δ*)    = {eps_opt:.4f}  (analítico: 0.2500)")
print(f"  M(δ*)    = {M_opt:.4f}  (analítico: 4.0000)")
print(f"  C(δ*)    = {C_opt:.4f}  (analítico: 4.0000)")

# ── Tabela 9 ─────────────────────────────────────────────────────────────────
print(f"\n--- Tabela 9: 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}")

#  Figura 4: Raio de convergência ε(δ), M(δ), C(δ) 
fig4, ax4 = plt.subplots(figsize=(7, 5))
ax4.plot(delta_curva, eps_curva, label=r"$\varepsilon(\delta)$",  color='steelblue', lw=2)
ax4.plot(delta_curva, M_curva,   label=r"$M(\delta)$",            color='tomato',    lw=2)
ax4.plot(delta_curva, C_curva,   label=r"$C(\delta)$",            color='seagreen',  lw=2)
ax4.axvline(d_opt,   color='k',         linestyle='--', lw=1.2, label=rf"$\delta^*={d_opt:.2f}$")
ax4.axhline(eps_opt, color='steelblue', linestyle=':',  lw=1.2, label=rf"$\varepsilon^*={eps_opt:.2f}$")
ax4.set_xlabel(r"$\delta$"); ax4.set_ylabel("")
ax4.set_ylim(0, 8); ax4.set_xlim(0, DELTA_MAX_CURVA)
ax4.set_title(r"Raio de convergência $\varepsilon(\delta)$, $M(\delta)$, $C(\delta)$")
ax4.legend(); ax4.grid(True, alpha=0.4)
plt.tight_layout()

plt.show()

#  Figura 5: 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]

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

plt.show()

#  Figura 6: N_max teórico vs cruzamento numérico 
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))
    mask = np.abs(t) <= 0.9
    cruzou_tp = next((n for n in range(N+1)
                      if np.max(np.abs(y_tp[n, mask] - x_ref[mask])) < tol), None)
    cruzou_rt = next((n for n in range(N+1)
                      if np.max(np.abs(y_rt[n, mask] - x_ref[mask])) < 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 (|t|≤0.9)':>18} {'n_rt (|t|≤0.9)':>18}")
print("-" * 58)
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:>18} {str(nrt)+flag_rt:>18}")
print("  † não atingiu a tolerância em N iterados")

fig6, ax6 = 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]
ax6.bar(x_pos - w, n_teo,    width=w, label=r"$N_{\max}$ teórico ($\delta^*$)", color='steelblue')
ax6.bar(x_pos,     n_num_tp, width=w, label="trapézio numérico",               color='tomato')
ax6.bar(x_pos + w, n_num_rt, width=w, label="retângulo numérico",              color='seagreen')
ax6.set_xticks(x_pos); ax6.set_xticklabels(labels)
ax6.set_xlabel(r"tolerância $\tau$"); ax6.set_ylabel("iterados necessários $n$")
ax6.set_title(r"(domínio $|t|\leq 0.9$; $\varepsilon_{\rm num}=0.99 > \varepsilon^*=0.25$)")
ax6.legend(); ax6.grid(True, axis='y', alpha=0.4)
plt.tight_layout()

plt.show()