# -*- coding: utf-8 -*-
"""
Created on Wed Oct  9 19:39:28 2024

@author: arnau
"""

# -*- coding: utf-8 -*-
"""
Created on Fri Oct  4 13:43:18 2024

@author: arnau
"""

import numpy as np
import matplotlib.pyplot as plt


#%% Gompertz

# constantes 

a = 0.036
kappa = 760
y0 = 16
T= 200

def solution(x):
    '''solution exacte de l'équation de Gompertz'''
    c = np.log(y0/kappa)
    return kappa*np.exp(c*np.exp(-a*x))

def Euler(y0,T,n):
    '''
    Entrées : [0,T] les extremités de l'intervalle de définition d'une solution
              y0 une condition initiale
              n le nombre de subdivision
    Sortie : une approximation de la solution de y'=F(y)
    sur [a,b] telle que y(a)=y0 avec Euler
    '''
    h=T/n  # pas 
    Y=[y0]
    for _ in range(1,n+1):
        yold=Y[-1]
        ynew = yold + h*a*yold*np.log(kappa/yold) 
        Y.append(ynew)# liste des approximation de y(a+k(b-a)/n)
    return Y

def Heun(y0,T,n):
    '''
    Entrées : [0,T] les extremités de l'intervalle de définition d'une solution
              y0 une condition initiale
              n le nombre de subdivision
    Sortie : une approximation de la solution de y'=F(y)
    sur [a,b] telle que y(a)=y0 avec Heun
    '''
    h=T/n  # pas 
    Y=[y0]
    for _ in range(1,n+1):
        ytemp = Y[-1]+h*a*Y[-1]*np.log(kappa/Y[-1])
        ynew = Y[-1] + h/2*(a*Y[-1]*np.log(kappa/Y[-1])+a*ytemp*np.log(kappa/ytemp) )
        Y.append(ynew)# liste des approximation de y(a+k(b-a)/n)
    return Y


#%% Tracé des solutions approchées et exactes

N = [1,5,10,100]
for n in N:
    temps = [k*T/n for k in range(n+1)]
    Y = Heun(y0,T,n)
    plt.plot(temps,Y, label = 'Heun avec n='+str(n))
    
plt.plot(np.linspace(0,T,1000), solution(np.linspace(0,T,1000)),':',label = 'solution exacte')
plt.legend()

#%% Croissance des rats musqués

Age =[ 0,21,29,35,44,50,55,62,69,76,83,90,97,104,110,117,124,132,138,146,152,180,187,194,201]  
Masse = [16,116,175,264,352,416,447,503,540,573,603,646,684,688,695,712,739,728,747,733,738,763,757,765,767 ]
plt.plot(np.linspace(0,T,1000), solution(np.linspace(0,T,1000)), 'blue',label = 'solution exacte')
plt.plot(Age,Masse, 'rx',label='Croissance du rat musqué')
plt.legend()
plt.show()


#%% Erreur

def Erreur(y0,T,n):
    Y1 = Euler(y0,T,n) # approximation par Euler
    Y2 = Heun(y0,T,n)  # approximation par Heun
    E1 = [np.abs(solution(k*T/n)-Y1[k]) for k in range(n+1)] # erreur pour Euler
    E2 = [np.abs(solution(k*T/n)-Y2[k]) for k in range(n+1)] # erreur pour Heun
    return max(E1),max(E2)

#%% erreur à T constant et pas qui tend vers $0$

Nb = range(1,100)  # nombre de subdivision
Erreur1= [Erreur(y0,T,k)[0] for k in Nb]
Erreur2= [Erreur(y0,T,k)[1] for k in Nb]
plt.plot(Nb,Erreur1, color='red', label='Erreur de la méthode Euler')  # erreur de la méthode d'Euler
plt.plot(Nb,Erreur2, color='blue',label='Erreur de la méthode de Heun') # erreur de la méthode de Heun
plt.title('Erreur lorsque le pas tend vers 0')
plt.xlabel('Nb de subdivision n')
plt.ylabel('Erreur')
plt.legend()

#%% erreur à pas constant avec T qui grandit

intervalle = range(1,50)
Erreur1= [Erreur(y0,k,100)[0] for k in intervalle]
Erreur2= [Erreur(y0,k,100)[1] for k in intervalle]
plt.plot(intervalle,Erreur1, color='red',label='Erreur de la méthode Euler')  # erreur de la méthode d'Euler
plt.plot(intervalle,Erreur2, color='blue',label='Erreur de la méthode de Heun') # erreur de la méthode de Heun
plt.title('Erreur à pas constant')
plt.xlabel('Longeur T de intervalle')
plt.ylabel('Erreur')
plt.legend()