"""
Thème 5 : Méthode d'Euler
"""

import numpy as np
import matplotlib.pyplot as plt


## Remarque préliminaire

# Remarque préliminaire : (E0) est une EDL1 à coefficients constants et homogène
#                          on sait la résoudre.
                          # Solutions : elle est unique puisqu'il y a une
                          # condition initiale ; c'est ys(t) = exp(t).

# ICI ON VA ETUDIER L'EQUATION DIFF. COMMME SI ON NE CONNAISSAIT PAS LA
# FONCTION EXP : LE BUT EST DE COMPRENDRE LA METHODE D'EULER A TRAVERS CET EXPLE.

## Q1 a)

# D'après (E0) y's(0)=ys(0) puisque ys vérifie l'éq. diff. Donc y's(0)=1

## Q1 -c)

# par déf, si f est solution, d'après (E0), on a f'(t)=f(t). Donc, on a la pente de la tangente
# à la courbe de f au point (t,f(t))


## Q1 -d)

# Une flèche F  sur le dessin indique que si une solution de y'=y passe par le point A d'où part cette fléche,
# alors la courbe de cette solution a pour tangente en A la flèche F


## Q2

# Etant donné le champ de vecteurs associé à une ED donnée, trouver une solution, c'est déterminer
# une courbe dont les tangentes sont prescrites, partant d'un point initial sur la courbe.

##Q3 Question classique d'oral ou d'écrit

# RAPPEL : y(tk) = valeur de la solution de (S), à la dae tk
#          yk    = valeur obtenue en tk en suivant le champ de vecteurs depuis le point initial de la courbe.
#          On a l'intuition que ce yk devrait donner une bonne approximation de y(tk).

# D'après l'intution on devrait avoir yk ~ y(tk) (~ signifie ici : "égale à peu près")
# D'après (S) y'(tk) = F(tk, y(tk)) donc  on a  y'(tk) ~  F(tk, yk) (*)

# Par aileurs, par défintion du taux d'accroissement et de la dérivée :
# y'(tk) ~ tx d'accroissement de y entre les instants tk et tk+1 si h = tk+1 - tk est petit
# donc on peut remplacer  dans (*),  y'(tk) par (yk+1 - yk)/ h, ce qui condiuit à :

#              (yk+1 - yk)/ h  ~  F(tk, yk) (**)

# il est donc raisonnable d'étudier la suite récurrente définie par :
#                          (yk+1 - yk)/ h  = F(tk, yk),
# dont les termes consécutifs devraient donner un nuage de points qui colle à la courbe de la solution y.








## RAPPEL : Programmation d'une suite récurrente
# exemple : un+1=f(un) et f(x)=1+x**2

def f(x):
    return 1+x**2

def suite(n):
    """
    renvoie la liste [u0,...,un]
    """
    u0 = 1     # mettons que u0 = 1
    Lu=[u0]    # liste des termes de la suite
    for k in range(n):
        x = Lu[-1]      # j'extraie le dernier terme de Lu
        y = f(x)        # je calcule son image par f
        Lu.append(y)     # je la mets en bout de liste
    return Lu



## Q5 Euler = programmation d'une suite récurrente

def EulerG(F,y0,a,b,n):
    h = (b-a)/n
    # intialisation de la liste des instants :
    Lt =[a] # a=t0
    # initialisation de la liste des valeurs approchées de y en tk
    Ly =[y0]

    tk,yk =Lt[-1],Ly[-1]  #  initialisation des variables tk et yk

    # Boucle de la suite récurrente  définie par (S') en page 18
    for k in range(n):
        yk += h*F(tk,yk)  # RAPPEL : x += a signifie x = x+a
        tk += h
        Lt.append(tk)
        Ly.append(yk)
    return Lt,Ly

## Q6

#a) (E0) est une EDL1 à coeff constants + condition initiale : elle admet une unique solution ye d'après le cours,
# elle est donnée par ye(t) = exp(t)

#b) on a vu qu'on obtient yk+1 = (1+h)yk puisque dans le cas
# de l'équation (E0), F(t,y) = y

#c) La suite (yk) est géométrique de raison r =(1+h) et
# de premier terme y0 = ye(0)= 1, donc yk = (1+h)**k (lire  : 1+h puissance k)


## Q7 Classique concours. Voir la solution rédigée sur papier
# Cette démo a été faite en classe (chap. 3bis)

# yN  : valeur approchée de ye(t_N)  (t indice N)

# Or, t_N = Nh, donc t_N = T

# Finalement : yN devrait être une valeur approchée de exp(T)

# preuve  : ln(yN)= N x ln (1 + T/N)  car h = T/N et yk = (1+h)^k

# par produit d'équivalents et équivalent usuel :

#          ln (yN) ~ T donc par composition de LIMITES yN -> exp(T)

# ainsi yN est bien une valeur approchée de ye(T) pour N grand.







## Q8a

def FE0(t,y):
    """
    fonction (t,y) --> y qui correspond à
    l'équa. diff (E0) y'=y
    """
    return y

# j'utilise Q5 puisque j'y ai programmé Euler en général
F,y0,a,b,n = FE0,1,0,3,30
Lt,Ly = EulerG(F,y0,a,b,n)
# solution exacte :
Lye = np.exp(Lt)   # j'utilise la vectorisation des fonctions numpy

plt.plot(Lt,Ly,'.--',label='solution approchée', )
plt.plot(Lt,Lye, label='solution exacte')
plt.grid('on')
plt.legend(loc='best')
plt.show()
# pour ceux qui veulent sauvegarder leur graphique :
plt.savefig('TP5-Q8a.pdf')

##Q9-10 Modèle logistique

# Même procédé qu'en Q8, puisqu'on a programmé EulerG :

def FL(t,y):
    """
    (t,y)--> y(1-y/10)
    """
    return y*(1-y/10)

def yL(t):
    return 10/(1+19*np.exp(-t))



F,y0,a,b,n = FL,0.5,0,5,100
# solution approchée
Lt,Ly = EulerG(F,y0,a,b,n)
# solution exacte :
Lyex = [yL(t) for t in Lt ]  #liste en compréhension

plt.plot(Lt,Ly,'.--',label='solution approchée', )
plt.plot(Lt,Lyex, label='solution exacte')
plt.grid('on')
plt.legend(loc='best')
plt.show()


##Q11 exemple non autonome


def FG(t,y):
    """
    (t,y)--> -2ty
    """
    return -2*t*y

# Remarque : (G) est  une EDL1 homogène, donc on
# sait calculer la solution exacte.


F,y0,a,b,n = FG,1,0,5,100
# solution approchée
Lt,Ly = EulerG(F,y0,a,b,n)

plt.plot(Lt,Ly,'.--',label='solution approchée', )
plt.grid('on')
plt.show()


## Q12

# Si je note ak la concentration en réactif A à la date tk, puis
# bk et ck ce que vous pensez avec des notations évidentes,
# le schéma d'Euler qu'on obtient est :


# ak+1 = ak -h x k1 x ak
# bk+1 = bk +h(k1 x ak -k2 x bk)
# ck+1 = ck +h x k2 x bk


## Q13
# La seule difficulté consiste à respecter une mise
# à jour synchrone des termes des suites (ak), (bk), (ck)
# pour garantir ce non-décalage temporel, FAIRE UNE AFFECTATION
# A LA VOLEE
a = 1   # mol/L
k1,k2 = 1,1
a0,b0,c0 = a,0,0
La,Lb,Lc,Lt =[a0],[b0],[c0],[0]

T = 6 # en secondes
n = 60
h = T/n # pas

for k in range(n):
    ak,bk,ck,tk = La[-1],Lb[-1],Lc[-1],Lt[-1]
    ak,bk,ck =  ak-h*k1*ak, bk+h*(k1*ak - k2*bk),  ck+h*k2*bk
    tk+=h
    La.append(ak)
    Lb.append(bk)
    Lc.append(ck)
    Lt.append(tk)

plt.plot(Lt,La,label='[A]')
plt.plot(Lt,Lb,label='[B]')
plt.plot(Lt,Lc,label='[C]')
plt.grid('on')
plt.legend(loc='best')
plt.show()