1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
| import matplotlib.pyplot as plt
from math import*
##Variables
x0 = 0.5 #valeur initiale de x
c0 = 0.5 #valeur initiale de c
r = 1.0 #coefficient d'erosion
rho = 0.1 #coefficient d'impatience
U = lambda x: -3 * exp(-0.38 *x) #fonction de contentement
U_der = lambda x: 3 * 0.38 * exp(-0.38 * x) #dérivée de U
D = lambda c: exp(c) - 2*c #fonction de sacrifice
D_der = lambda c: exp(c) - 2 #dérivée de D
phi = lambda l: log(l + 2) if l >= D_der(0) else 0 #phi, inverse de D'
t_max = 5.0 #temps maximal jusqu'auquel on regarde
n = 1000 #nombre d'itérations
##Differential equation system
h = t_max/n #le pas temporel
X = [x0] #la liste des x(t)
L = [D_der(c0)] #la liste des l(t), l(t) désignant lambda(t) dans le sujet
for i in range(n) :
x = X[i] #on prend la valeur actuelle de x(t)
l = L[i] #et celle de l(t)
x_der = -r * x + phi(l) #on calcule x'(t)
l_der = -(r + rho) * l + U_der(x) #et l'(t)
X.append(x + h * x_der) #et on prend x(t+h) = x(t) + h * x'(t) (méthode d'euler)
L.append(l + h * l_der) #pareil pour l(t)
C = [phi(l) for l in L] #on prend la liste des efforts qui valent phi appliquées aux l(t)
plt.figure() #et on trace
plt.plot([h * i for i in range(n + 1)], X)
plt.title('Evolution du bien-être du couple')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.figure()
plt.plot([h * i for i in range(n + 1)], C)
plt.title("Evolution de l'effort")
plt.xlabel('t')
plt.ylabel('c(t)')
plt.show() |
Partager