from math import *
import matplotlib.pyplot as plt

from scipy.integrate import odeint

## Q1
# RFD :l'accélération intervient souvent, avec des forces
# ou des moments qui dépendent de la position

## Q2
# Y est un vecteur 2D maintenant
def Fvolterra(Y,t):
    (x,y)=Y
    dx=x*(2/3-4/3*y)
    dy=y*(x-1)
    return (dx,dy) # ou éventuellement [dx,dy]

N=500
tMax=30
tList=[i*tMax/(N-1) for i in range(N)] # on peut diviser par N...

## Q3
# renvoie la liste des [x,y], pour les 3 dates demandées

## Q4
# on peut définir une fonction avec un argument optionnel
def resoutVolterra(x0,y0,solveur=odeint):
    # N ettList sont connues (var globales)
    # on peut coder ainsi pour éviter les "append" : un peu moins coûteux
    xList=[0]*N
    yList=[0]*N # donne [0,0,...,0] N fois

    # 3 arguments, les CI sont un vecteur 2D
    # solveur vaut PAR DEFAUT odeint
    toutLeBazar = solveur(Fvolterra,(x0,y0),tList)

    for i in range(N): # la longueur est connue
        xList[i],yList[i]=toutLeBazar[i] # on peut aussi coder ...[i][0]

    return (xList,yList)

## Q5
# Il faut bien sûr appeler la fonction codée
xList,yList=resoutVolterra(2,0.6)
# 2 courbes = 2 appels à plot, le show() à la fin
# on doit appeler plt.legend() pour afficher les labels
plt.plot(tList,xList,'b-',label='Proies')
plt.plot(tList,yList,'r-',label='Prédateurs')
plt.legend()
plt.show()

## Q6
# Les arguments obligatoires sont : nomDeFn, CI, listeDates
# on doit renvoyer la liste des [x,y]
def euler2D(F,Y0,tList):
    # pas forcément le même tList que la var globale :
    # sa longeur n'est pas connue
    N=len(tList)
    res=[Y0]
    Yk=Y0 # variable de boucle
    for k in range(N-1):
        xk,yk=Yk
        tk=tList[k]
        tkp1=tList[k+1]
        dt=tkp1-tk
        dx,dy=F(Yk,tk)
        # on modifie la var de boucle
        Yk=[xk+dx*dt, yk+dy*dt]
        # et on la stocke surtout !!
        res.append(Yk)
    return res

## Q7 : remplacer simplement odeint par euler2D ligne 33
## ou appeler ligne 43 : xList,yList=resoutVolterra(2,0.6,solveur=euler2D)