#!/usr/bin/python3
#
# Un deuxième pour la méthode d'Euler: un objet ralenti par le frottement

# Constantes physiques:
alpha = 1.  
beta = 2.
v0 = 3.

# Constantes numériques: 
T = 8.     # intervalle de temps
N = 1000   # nombre de pas à calculer
h = T / N  # pas d'incrément

# La fonction du membre de droite de l'EDO
def f(v, t):
    return -alpha * v - beta * v**2
        
# Importer les bibliothèques        
from euler import euler
import numpy as np
import matplotlib.pyplot as plt

# Calculer la solution
t = np.linspace(0, T, N+1)
solution = euler(f, v0, h, N)

# Affichage
plt.rc('text', usetex=True)
plt.plot(t, solution, label=r"solution numérique") # tracer la solution
plt.plot(t, v0 / (1 + beta * v0 * t), 'g--',label="approximation de grande vitesse")
tgtr1 = t[125:] # t >= 1
v1 = solution[125]
plt.plot(tgtr1, v1 * np.exp(-alpha*(tgtr1-1)), 'r--', label="approx. de petite vitesse à partir de t=1")
plt.xlabel("$t$", fontsize=20)
plt.ylabel("$v(t)$", fontsize=20)
plt.legend(fontsize=20)
plt.figure() # Nouvelle figure pour tracer l'erreur numérique
plt.plot(t, solution - alpha * v0 / (np.exp(alpha * t) * (alpha + beta * v0) - beta * v0))
plt.xlabel(r"$t$",fontsize=20)
plt.ylabel(r"$v_{{\rm num}}(t) - v(t)$",fontsize=20)
plt.figure() # erreur numérique relative
plt.plot(t, (np.abs(solution - alpha * v0 / (np.exp(alpha * t) * (alpha + beta * v0) - beta * v0))) \
              / (alpha * v0 / (np.exp(alpha * t) * (alpha + beta * v0) - beta * v0)))
plt.xlabel(r"$t$",fontsize=20)
plt.ylabel(r"$|v_{{\rm num}}(t) - v(t)| / |v(t)|$",fontsize=20)

plt.show()


