# -*- coding: utf-8 -*-
''' TP 01 proposé par Magali
ce sont mes réponses
'''

import numpy as np
import matplotlib.pyplot as plt


#%% premier exo
K=70
U0=1.8

def f0(t,y):
    return -K*y



def g0(t):
    return U0*np.exp(-K*t)


def euler(f,y0,a,b,n):
    """méthode d'Euler
        f fonction défininssant l' équation diff
         y0 valeur initiale en a
         b borne de l'intervalle de calcul 
         n nombre de pas
    """
    pas= (b-a)/n
    resultat=[y0]
    y=y0
    t=a
    for i in range(1,n+1):
        y=f(t,y)*pas +y
        resultat.append(y)
        t=t+pas
    return resultat         



n=30
Y=euler(f0, U0, 0, 0.5, n)
T=np.linspace(0,0.5,n+1)
plt.plot(T,Y)
plt.plot(T,g0(T))
plt.show()

#%%


for n,c in [(20,"red"),(30,"goldenrod"),(60,"teal"),(160,"darkorchid")]:
    Y=euler(f0, U0, 0, 0.5, n)
    T=np.linspace(0,0.5,n+1)
    plt.plot(T,Y)
    plt.plot(T,g0(T),color=c)
plt.show()
#%% deuxième exercice

def f1(t,y):
    return y+2*t

def g1(t):
    return 3*np.exp(t)-2*t-2



n=40
Y=euler(f1, 1, 0, 4, n)
T=np.linspace(0,4,n+1)
plt.plot(T,Y)
plt.plot(T,g1(T))
plt.show()
#%% mesure de l'erreur max

def ecartMax(L1,L2):
    if len(L1)!=len(L2) or L1==[]:
        raise ValueError
    else:
        m=abs(L1[0]-L2[0])
        for i in range(len(L1)):
            m=max(m,abs(L1[i]-L2[i]) )
    return m

#%%

for n in [10,20,50,100,300,4000,5000]:
    Y=euler(f1, 1, 0, 4, n)
    T=np.linspace(0,4,n+1)
    print(ecartMax(Y,g1(T)))
    
#%% Méthode de Heun
def heun(f,y0,a,b,n):
    """méthode d'Heun
        f fonction défininssant l' équation diff
         y0 valeur initiale en a
         b borne de l'intervalle de calcul 
         n nombre de pas
    """
    pas= (b-a)/n
    resultat=[y0]
    y=y0
    t=a
    for i in range(1,n+1):
        y1=f(t,y)*pas +y
        y=((f(t,y)+f(t+pas,y1))/2)*pas+y
        resultat.append(y)
        t=t+pas
    return resultat         
#%%
for n in [10,20,50,100,300,4000,5000]:
    Y=euler(f1, 1, 0, 4, n)
    Z=heun(f1, 1, 0, 4, n)
    T=np.linspace(0,4,n+1)
    print(ecartMax(Y,g1(T)),"     " ,ecartMax(Z,g1(T)))