from math import *
import matplotlib.pyplot as plt

from scipy.integrate import odeint

## Méthode d'Euler pour les équations différentielles d'ordre 1
## On travaille dans tout le TD sur l'équadiff

##    τ y' + y = ...

## Dans le TD, on fixe la constante de temps τ=1
## et on prendra toujours la valeur initiale y0=1

## Liste des dates
# nombre de points
N=80
## Q5
## Créer la variable tList contenant N points de 0 à 5 inclus, régulièrement espacés
pas=5/(N-1)
tList=[i*pas for i in range(N)]

## Q6
## Compléter : c'est la fonction F fondamentale de l'équadiff qui calcule
## la dérivée y'; ici l'équadiff est avec excitation.
def Fequadiff(y,t):
    return -y+cos(8*t)

## Q7
## Compléter pour retourner la liste [yk] nommée yList, impérativement de même longueur
## que tList (sinon, une erreur surviendra lors du tracé).
## (c'est une bonne idée de le vérifier après l'appel, grâce à la commande len).
# Construit la liste des estimations des yk : fonction sans argumement car la valeur
# initiale est fixée à 1, et car tList et la fonction de l'équadiff sont des variables
# globales.

def resoutEuler():
    yk=1
    yList=[yk]
    for i in range(N-1): # si N, tList[i] déclenchera une erreur à la fin
        tk=tList[i]
        tkp1=tList[i+1]
        yk=yk+Fequadiff(yk,tk)*(tkp1-tk)
        yList.append(yk)
    return yList

## Q8
## Expliquer les deux premières lignes du code : voir l'aide de odeint
## dans le document d'accompagnement.
def resoutOde():
    res=odeint(Fequadiff,1,tList)
    yList=[r[0] for r in res]
    return yList

yList=resoutEuler()
plt.plot(tList,yList,'b-',label='Euler')
yList=resoutOde()
plt.plot(tList,yList,'g-',label='ODE(rk4)')
## Q9
## solution de l'équadiff avec excitation
def fonctionVraie(t):
    return (1-cos(atan(8))/sqrt(65))*exp(-t)+1/sqrt(65)*cos(8*t-atan(8))

yList=[fonctionVraie(t) for t in tList]
plt.plot(tList,yList,'r.',label='exact')
plt.legend()
plt.show()