#### Equations différentielles - Méthode d'Euler (ou Point Milieu)
## Armel MARTIN - Septembre 2021
## distribué sous licence CC-by-nc-sa
import numpy as np
import matplotlib.pyplot as plt
#from scipy import integrate


# second membre du problème de cauchy 
def cauchy_1d(y,t):
        '''Calcule le second membre du problème de Cauchy 1D: y' = f(y,t)'''
        return -y+1 # équation y' = -y+1
#        return -y+t # équation y' = -y+t
#        return -y+np.cos(2*np.pi*t) # équation y' = -y+t

def euler_1d(f,y0,t0,t1,N):
        '''Résoud le problème de Cauchy 1D (ED du 1er Ordre) y' = f(y,t):
        _avec f donnée en argument,
        _sur l'intervalle [t0,t1],
        _avec pour condition initiale y(t0) = y0,
        _en N pas de temps.
        Renvoie: les vecteurs des instants et des valeurs de y à ces instants.'''
        h = (t1-t0)/N # définition du pas de temps
        ## méthode avec des listes
        t = t0
        y = np.copy(y0) # préférable de copier proprement pour que y0 ne soit pas
                     # modifiée si on s'en ressert un peu plus loin dans la fonction
        list_t = [t0]
        list_y = [y0]
        for i in range(N):
                y = y + h*f(y,t) # Pas d'Euler
                #y = y + h*f( y + h/2*f(y,t) , t + h/2) # Pas du point milieu
                list_y = list_y + [y]
                t = t + h
                list_t = list_t + [t]
        return np.array(list_t),np.array(list_y) # on renvoie des array pour plot ensuite

## méthode avec des "array" (classe ndarray):
## #attention cette méthode nécessite de connaître a priori la taille des tableaux qu'on va créer,
## ce qui est le cas ici car on utilise un pas de temps fixe.
##	t = arange(t0,t1+h,h) # N+1 points
##	y = zeros(N+1)
##	y[0] = y0
##	for i in range(N):
##		t[i+1] = t[i] + h
##		y[i+1] = y[i] + h*f(y[i],t[i]) # pas d'Euler
##	return t,y

# Condition initiale et domaine temporel:
y0,t0,t1,N = 0,0,8.,500
# Graphe
plt.figure(1,figsize=[4,3])
sol_euler = euler_1d(cauchy_1d,y0,t0,t1,N) # résolution numérique
yy = sol_euler[1] # solution numérique Euler
tt = sol_euler[0] # axe des temps
plt.plot(tt,yy,label='N='+str(N))
tt = np.linspace(t0,t1,N)
yy_ex = 1+(y0-1)*np.exp(-tt) # Solution exacte: y' = -y+1
#yy_ex = tt-1+(y0+1)*np.exp(-tt) # Solution exacte: y' = -y+t
#C = np.cos(np.arctan(2*np.pi))/np.sqrt(1+4*np.pi**2)
#yy_ex = np.cos(2*np.pi*tt-np.arctan(2*np.pi))/np.sqrt(1+4*np.pi**2) + \
#        (y0-C)*np.exp(-tt) # Solution exacte: y' = -y+cos(2*pi*t)
plt.plot(tt,yy_ex,color='black',label='exacte')
plt.xlabel('temps (unité tau)')
plt.ylabel('y')
#plt.xlim((0,14))
#plt.ylim((-2.5,2.5))
plt.title('dy/dt = -y+1')
#plt.title('dy/dt = -y+t')
#plt.title('dy/dt = -y+cos(2*pi*t)')
plt.legend(loc='upper right')
plt.grid()
plt.tight_layout()
plt.savefig('cn3_2.png')
#plt.savefig('cn3_3.png')
#plt.savefig('cn3_4.png')
plt.show()