# -*- coding: utf-8 -*-
"""
@author: E. Lins
"""
import numpy as np
import matplotlib.pyplot as plt

# Funcao que retorna a solucao exata
def f_exata(t):        
    return ( (1+t)**2 - 0.5*np.exp(t))

# Funcao que define o EDO
def f(t, w):
    return (w - t**2 + 1)

# Metodo Runge Kutta de 4a ordem
def RungeKutta4(t0, y0, t_vals, h, N):
    y_vals = np.zeros(N+1)
    y_vals[0] = y0
    
    for i in range(N):
        t = t_vals[i]
        y = y_vals[i]
        
        k1 = h * f(t, y)
        k2 = h * f(t + h/2, y + k1/2)
        k3 = h * f(t + h/2, y + k2/2)
        k4 = h * f(t + h, y + k3)
        
        y_vals[i+1] = y + (k1 + 2*k2 + 2*k3 + k4) / 6
    
    return y_vals

# Metodo Adams-Bashforth de 4a ordem
def AdamsBashforth4(t0, y0, t_vals, h, N):
    y_vals = np.zeros(N+1)
    y_vals[0] = y0

    # Calcula os valores iniciais usando o RK4
    y_start = RungeKutta4(t0, y0, t_vals[:4], h, 3)
    y_vals[:4] = y_start[:4]

    # Coeficientes para o método de Adams-Bashforth de 4a ordem
    a1, a2, a3, a4 = 55/24, -59/24, 37/24, -9/24

    for i in range(3, N):
        t = t_vals[i]
        y = y_vals[i]        
        
        # Previsão usando Adams-Bashforth de 4ª ordem
        y_vals[i+1] = y_vals[i] + h * (a1*f(t, y) + a2*f(t_vals[i-1], y_vals[i-1]) +
                                       a3*f(t_vals[i-2], y_vals[i-2]) + a4*f(t_vals[i-3], y_vals[i-3]))
    
    return y_vals

# Entradas: intervalo e condições iniciais
a = 0  # Limite inferior do intervalo
b = 2  # Limite superior do intervalo
N = 10  # Número de subintervalos (passos)
alpha = 0.5  # Valor inicial da solucao no ponto t=a

h = (b-a)/N  # Tamanho do passo
t = np.linspace(a, b, N+1)  # Gera os pontos t no intervalo [a, b]

# Solucao com Adams-Bashforth 4ª ordem
ab4 = AdamsBashforth4(a, alpha, t, h, N)

# Calcula a solucao exata para comparacao
te = np.linspace(a, b, 10*N)  # Gera os pontos t com maior refinamento para solucao exata
fe = f_exata(te)  # Solucao exata

plt.plot(te, fe, label="Exata")  
plt.plot(t, ab4, 'r+:', label="Adams-Bashforth 4a ordem")  

# Configurações do gráfico
plt.xlabel("t")  # Rótulo do eixo t
plt.ylabel("y")  # Rótulo do eixo y
plt.grid()  # Adiciona grade 
plt.legend()  # Adiciona legenda