#########################################
######## TP informatique bio-sép ########
########        2025-2026        ########
#########################################


#########################################
## Thème 5 : Equations différentielles ##
#########################################

import numpy as np
import matplotlib.pyplot as plt

## Q1 : interprétation d'une ED1 avec condition initiale

#Q1a : y'(0) = y(0) = 1

#Q1b : la tangente à la courbe représentative de y
#       a pour pente 1 et passe par le point (0,1).
#      Son équation est donc : y = x+1

#Q1c : la tangente au point (t,f(t)) a pour pente f'(t) = f(t)

#Q1d : le champ de vecteurs représente en chaque point
#      la tangente à la courbe d'une solution de l'ED

## Q2 : Reformulation du problème

# Etant donné un champ de vecteurs associé à une équation différentielle,
# trouver une solution c'est déterminer une courbe dont les tangentes
# sont prescrites, partant d'un point initial sur la représentation
# du champ de vecteurs.

## Q3 : Justification de la méthode d'Euler

# Si y est dérivable en tk, on a le DL d'ordre 1 :
# y(tk+h) = y(tk) + h*y'(tk) + o(h)
#
# En posant t_{k+1} = tk + h :
# y(t_{k+1}) ~ y(tk) + h*y'(tk)
#            ~ y(tk) + h*F(tk,y(tk))

## Q4 : Calcul des premiers termes par la méthode d'Euler

# I = [0,3] et n=6 donc on discrétise l'intervalle I
# à l'aide des valeurs :
# t0 = 0, t1 = 0.5, t2 = 1, t3 = 1.5 ... t6 = 3
# c'est-à-dire : tk = a + k*(b-a)/n pour tout k dans [[0,6]]
# en posant : I = [a,b] donc a = 0, b = 3
#
# L'approximation précédente donne ici :
# y(t_{k+1}) ~ y(tk) + 0.5*y(tk)
#            ~ 1.5*y(tk)
#
# d'où le tableau :
#
# k = | tk = | yk = | y'(tk)= |
#   0 |  0   |  1   |   1     |
#   1 |  0.5 |  1.5 |   1.5   |
#   2 |  1   |  2.25|   2.25  |
#   3 |  1.5 |  3.37|   3.37  |
#   4 |  2   |  5.06|   5.06  |
#   5 |  2.5 |  7.59|   7.59  |
#   6 |  3   | 11.39|  11.39  |

## Q5 : Méthode d'Euler (à gauche) pour une ED d'ordre 1

def EulerG(F,y0,a,b,n) :
    T = [a]
    Y = [y0]
    h = (b-a)/n
    for _ in range(n) :
        y0 = y0 + h*F(a,y0)
        Y.append(y0)
        a = a + h
        T.append(a)
    return T,Y

def F(t,y) :
    return y

## Q6/Q7/Q8 : Vérification sur un exemple

# Q6abc. La solution exacte est : ye(t) = exp(t)

# Q7. Sur [0,T], on a h = (T-0)/N = T/N
#     et la discrétisation de la méthode d'Euler
#     conduit à : y_{k+1} = yk + h*yk
#                         = (1+h)*yk
#     La suite (yk) est donc géométrique de raison (1+h).
#     On a donc : pour tout k, yk = y0*(1+h)**k
#     Pour k = N, on obtient : yN = 1*(1+h)**N
#                                 = (1+h)**(T/h)
#     La question est donc de trouver la limite
#     quand h tend vers 0 de : (1+h)**(T/h)
#     On utilise la forme exponentielle :
#         (1+h)**(T/h) = exp( T/h * log(1+h) )
#     et un équivalent en 0 : log(1+h) ~ h
#     donc T/h * log(1+h) ~ T/h * h = T
#     donc T/h * log(1+h) tend vers T
#     et par composée de limites,
#        (1+h)**(T/h) tend vers exp(T).
#     CQFD.

# Q8. Comparaison graphique

for n in [6,30,300,3000] :
    T,Y = EulerG(F,1,0,3,n)
    l = 'n = '+str(n)
    plt.plot(T,Y,'-.',label = l)

X = np.linspace(0,3,1000)
Y2 = np.exp(X)
plt.plot(X,Y2,color = 'pink', label = "y = exp(t)")
plt.title("Méthode d'Euler pour l'équation de Malthus")
plt.legend()
plt.show()

## Q9/Q10 : Vérification sur le modèle logistique

def FLogis(t,y) :
    return y*(1-y/10)

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

for n in [6,30,300,3000] :
    T,Y = EulerG(FLogis,0.5,0,10,n)
    l = 'n = '+str(n)
    plt.plot(T,Y,'-.',label = l)

X = np.linspace(0,10,1000)
Y2 = solution(X)
plt.plot(X,Y2,color = 'pink', label = r"$y = \frac{10}{1+19exp(-t)}$")
plt.title("Méthode d'Euler pour le modèle logistique")
plt.legend()
plt.show()

## Q11 : Une équation différentielle non autonome

def F11(t,y) :
    return -2*t*y

T,Y = EulerG(F11,1,0,10,1000)

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

Yexact = [solution_exacte(t) for t in T]

plt.plot(T,Y,'-.',label = 'valeur approchée')
plt.plot(T,Yexact,color = 'pink', label = 'solution exacte')
plt.title("Exemple avec une ED non autonome")
plt.legend()
plt.show()

## Q12 : un système différentiel

def EulerSysteme(a,k1,k2,n) :
    '''resout le systeme differentiel sur [0,6] avec n points de calcul'''
    Lt = [0]
    La = [a] # concentration initiale en produit A
    Lb = [0] # concentration initiale en produit B
    Lc = [0] # concentration initiale en produit C
    h = 6/n  # résolution sur [0,6]
    for _ in range(n) :
        a0,b0,c0 = La[-1], Lb[-1], Lc[-1]
        a1 = a0 + h*(-k1*a0)
        b1 = b0 + h*(k1*a0-k2*b0)
        c1 = c0 + h*(k2*b0)
        La.append(a1)
        Lb.append(b1)
        Lc.append(c1)
        Lt.append(Lt[-1]+h)
    return Lt,La,Lb,Lc

## Q13 : tests avec plusieurs valeurs de k1 et k2

# k1,k2 = 1,1
# k1,k2 = 10,1
k1,k2 = 1,10
a = 1
n = 1000
T,A,B,C = EulerSysteme(a,k1,k2,1000)
plt.plot(T,A,label = 'produit A')
plt.plot(T,B,label = 'produit B')
plt.plot(T,C,label = 'produit C')
plt.legend()
plt.show()