changement de git

EDO.py

+#Etape 1 :  Écrire un solveur pour la méthode d'Euler explicite (méthode de Runge--Kutta à un étage), un pour la méthode de Runge--Kutta à deux étages et un pour la méthode de Runge--Kutta d'ordre 4.
+
+import numpy as np
+from functools import partial
+import matplotlib.pyplot as plt
+from mpl_toolkits import mplot3d
+
+# dyx/dt=σ(yy−yx),
+# dyy/dt=yx(ρ−yz)−yy,
+# dyzdt=yy*yx−βyz,
+
+def lorenz(y, t, sig, rho, beta):
+    return np.array([sig * (y[1]-y[0]), y[0] * (rho - y[2]) - y[1], y[1]*y[0] - beta * y[2]])
+
+def rk1(y, t, dt, F): # Euler
+    return y + dt * F(y,t)
+
+def k1(y, t, dt, F):
+    return dt * F(y, t)
+    
+def k2(y, t, dt, F):
+    return dt * F(y + (k1(y, 0, dt, F) / 2), t + (dt / 2)) 
+
+def k3(y, t, dt, F):
+    return dt * F(y + (k2(y, 0, dt, F) / 2), t + (dt / 2)) 
+
+def k4(y, t, dt, F):
+    return dt * F(y + (k3(y, 0, dt, F)), t + dt) 
+
+def rk4(y, t, dt, F):
+    return y + (1 / 6) * (k1(y, 0, dt, F) + 2 * (k2(y, 0, dt, F) + k3(y, 0, dt, F)) + k4(y, 0, dt, F))
+"""
+sig = 10
+beta = 8/3
+rho = 28
+dt = 0.001
+t_max = 1.5586522
+
+y1 = np.array([-0.9101673912,-1.922121396,18.18952097]) # état initial
+y2 = np.array([-0.9101673912,-1.922121396,18.18952097]) # état initial
+
+xRK1 = np.array([])
+yRK1 = np.array([])
+zRK1 = np.array([])
+
+xRK4 = np.array([])
+yRK4 = np.array([])
+zRK4 = np.array([])
+
+for i in range(int(t_max / dt)):
+    y1 = rk1(y1, 0, dt, partial(lorenz, sig=sig, beta=beta,rho=rho))
+    xRK1 = np.append(xRK1, y1[0])
+    yRK1 = np.append(yRK1, y1[1])
+    zRK1 = np.append(zRK1, y1[2])
+
+    y2 = rk4(y2, 0, dt, partial(lorenz, sig=sig, beta=beta,rho=rho))
+    xRK4 = np.append(xRK4, y2[0])
+    yRK4 = np.append(yRK4, y2[1])
+    zRK4 = np.append(zRK4, y2[2])
+
+ax = plt.axes(projection='3d')
+ax.plot3D(xRK1, yRK1, zRK1, 'blue')  
+ax.plot3D(xRK4, yRK4, zRK4, 'green')
+plt.show()
+"""
+
+#Etape 2 : Équation de Lorenz
+#Afin de valider les solveurs, appliquer les trois solveurs à l'équation de Lorenz et reproduire une figure approchant celle de l'énoncé (faire un graphique éventuellement tri-dimmensionnel de la trajectoire obtenue pour une condition initiale quelconque).
+
+# Etape 3 : Orbites périodiques
+"""L'attracteur de Lorenz contient des trajectoires périodiques sur
+lesquelles le système revient au point initial après un certain temps.
+Ces trajectoires sont dites "instables", comme toute autre trajectoire
+de ce système : il suffit qu'on dévie un tout petit peu de la
+trajectoire, et on finit par s'en éloigner complètement. Le point
+suivant se trouve (dans la limite de la précision numérique fournie) sur
+une trajectoire périodique de période $t_{max}=1.5586522$:
+$y_0=(-0.9101673912,-1.922121396,18.18952097)$. Prendre ce point comme
+valeur initiale, et faites évoluer le système avec un pas de temps donné
+(par exemple $\delta t = 0.001$). Combien de temps arrivez-vous à rester
+sur l'orbite périodique avec chacun des solveurs? Tracez les deux
+trajectoires sur une figure 3D pour les comparer.
+
++
+"""
+#Etape 4 : Modélisation d'épidémies
+"""Une fois les solveurs validés vous vous attaquerez à la tâche d'utiliser vos solveurs
+pour résoudre un vrai problème d'actualité.
+"""
+
+#Etape 5 : Le modèle SEIR
+
+def epidemie(y, t, R0, Tinf, Tinc):
+
+    return np.array([-(R0 / Tinf) * y[2] * (y[0] / N), (R0 / Tinf) * y[2] * (y[0] / N) - (1 / Tinc) * y[1],
+            (1 / Tinc) * y[1] - (1 / Tinf) * y[2], (1 / Tinf) * y[2]])
+
+
+y1 = np.array([9900, 100, 0, 0]) # SEIR
+
+R0 = 100
+Tinf = 15
+Tinc = 5
+N = 10000
+
+dt = 0.1
+t_max = 200
+t = 0
+t_list = np.array([])
+
+S = np.array([])
+E = np.array([])
+I = np.array([])
+R = np.array([])
+
+for i in range(int(t_max / dt)):
+    y1 = rk4(y1, 0, dt, partial(epidemie, R0=R0, Tinf=Tinf, Tinc=Tinc))
+    print(y1)
+    S = np.append(S, y1[0])
+    E = np.append(E, y1[1])
+    I = np.append(I, y1[2])
+    R = np.append(R, y1[3])
+
+    t += dt
+    t_list = np.append(t_list, t)
+
+    Nt = y1[0] + y1[1] + y1[2] + y1[3]
+
+    print("S : " + str(y1[0]))
+    print("E : " + str(y1[1]))
+    print("I : " + str(y1[2]))
+    print("R : " + str(y1[3]))
+    print(N)
+    print(Nt)
+
+plt.plot(t_list, S, 'b')
+plt.plot(t_list, E, 'r')
+plt.plot(t_list, I, 'g')
+plt.plot(t_list, R, 'purple')
+
+plt.legend(['S', 'E', 'I', 'R', 'swiss'])
+plt.show()
+
+#Etape 6 : cas reel
\ No newline at end of file