##Q1 euler def euler (n): h = (b - a) / n T=[a] t=a V=[0.]*n V[0] = 0 for i in range(n-1): V[i+1]=V[i]+h*F(V[i]) t += h T.append(t) #Q3 #plusieurs tactique import numpy as np pas = Tmax / (N - 1) N = int (Tmax / pas ) + 1 T = np . linspace(0 , Tmax , N) T = np . arange (0 , Tmax + pas , pas ) def Tableau ( tmax , h ) : return ( np . arange (0 , tmax + h , h ) ) T = Tableau (Tmax , pas ) #Q4 def f1(ti, yi) : # Cette fonction ne dépend pas directement de ti car elle est linéaire en yi return( np.array([yi[1], pow(omega0, 2) * (K - yi[0]) - 2 * m * omega0 * yi[1]]) ) #Q6 def EulerExplicite(Yini, h, Tmax, F) : '''renvoie le tableau de la solution de l'équation dY/dt = F(t, Y).''' T = Tableau(Tmax, h) SY = np.zeros([len(T), 2]) # tableau de type array à coefficients nuls SY[0] = Yini # première ligne du tableau renseignée for i in range(len(T) - 1) : # on s'arrête à l'avant dernier terme SY[i + 1] = SY[i] + h * F(T[i], SY[i]) return( SY ) ##pendule avec EUler def F(y,t): return np.array([y[1],-y[1]*ca-k*np.sin(y[0])]) T=np.linspace(0,tm,int(tm/h)+1) y0=np.array([theta0,v0]) Y=odeint(F,y0,T) pl.plot(T,np.array(Y)[:,0],label=str(theta0)) #Q6. ici on travaille avec le vecteur solution de l'équation différentielle Y'=F(Y) avec Y la matrice dont les lignes sont les instants t_i et les deux colonnes [theta,v] #Q7. pendule modifié import matplotlib.pyplot as pl import numpy as np from scipy.integrate import odeint xi = int(input("Entrer le facteur d'amortissement : ")) omega0 = int(input("Entrer la pulsation propre : ")) theta0 = int(input("Entrer la position initiale : ")) v0 = float(input("Entrer la vitesse angulaire initiale : ")) tm =int( input("Entrer la periode de simulation : ")) h =float( input("Entrer le pas d'echantillonnage : ")) k=omega0**2 ca=2*xi*omega0 def F(y,t): return np.array([y[1],-y[1]*ca-k*np.sin(y[0])]) # Graphe de x en fonction de t pl.subplot(2,1,1) c,d=-3,25 # tm=120 # h=1/500 T=np.linspace(0,tm,int(tm/h)+1) y0=np.array([theta0,v0]) Y=odeint(F,y0,T) pl.plot(T,np.array(Y)[:,0],label=str(theta0)) pl.grid() # Grille pl.title("Position") pl.xlabel("t") pl.ylabel("theta(t)") pl.axis([0,tm,c,d]) pl.legend(loc="lower right") #Q8. Portrait de phase pl.subplot(2,1,2) # y0=np.array([0,1]) a,b=-3,20 c,d=-2,3 y0=np.array([theta0,v0]) Y=odeint(F,y0,T) pl.plot(np.array(Y)[:,0],np.array(Y)[:,1],label=str(theta0)) pl.grid() # Grille pl.title("Portraits de phase") pl.xlabel("theta(t)") pl.ylabel("theta'(t)") pl.axis([a,b,c,d]) pl.legend(loc="lower right") pl.show() ## Autoroute A7 import matplotlib.pyplot as plt import numpy as np from math import pi ### carte def trace(q_exp, v_exp): nb_mesures = len(q_exp) c_exp = [0] * nb_mesures for i in range(nb_mesures): c_exp[i] = q_exp[i] / v_exp[i] plt.plot(c_exp, q_exp, 'o') plt.xlabel('Concentration (véhicule/km)') plt.ylabel('Débit (véhicule/h)') plt.show() q_exp=np.array([ 0., 2000, 4000, 6000, 7000, 6000, 7000, 6000, 7000, 6500]) v_exp=np.array([ 100, 100, 100, 100, 100, 100, 100, 80, 50, 40 ]) trace(q_exp,v_exp) ### question 1 def debit(v_max, c_max, C_ligne): debit = [] for concentration in C_ligne: vitesse = v_max * (1 - concentration / c_max) debit.append(concentration * vitesse) return debit def diagramme(v_max, c_max, C_ligne): q = debit(v_max, c_max, C_ligne) plt.plot(C_ligne, q, marker='o') plt.xlabel('Concentration (véhicule/m)') plt.ylabel('Débit (véhicule/s)') plt.plot((0.04, c_max/2, c_max/2), (c_max * v_max / 4, c_max * v_max / 4, 0.2),linestyle=":") plt.text(0, c_max * v_max / 4, r'$\frac{c_{max} . v_{max}}{4}$', verticalalignment='center', fontsize=15) plt.text(c_max/2, 0, r'$\frac{c_{max}}{2}$', horizontalalignment='center', fontsize=15) plt.show() return q v_max, c_max = 130*1000/3600, 200/1000 diagramme(v_max, c_max, np.linspace(0,c_max)) ### question 4 def C_depart(La, dx, c1, c2, d1, d2): C0 = [] for l in np.arange(0, La, dx): if l < d1 or l > d2: C0.append(c1) else: C0.append(c2) return C0 ### question 6 def resolution(C , dt , dx , c_max , v_max): for i in range(1,len(C)): Q = debit(v_max , c_max ,C[i-1]) Q.append(Q[0]) # ajoute le véhicule de droite du dernier for j in range(len(C[i])): C[i][j] = C[i-1][j] - dt * (Q[j+1] - Q[j]) / dx def resolution_arriere(C , dt , dx , c_max , v_max): for i in range(1,len(C)): Q = debit(v_max , c_max ,C[i-1]) Q.append(Q[0]) # ajoute le véhicule de droite du dernier for j in range(len(C[i])): C[i][j] = C[i-1][j] - dt * (Q[j] - Q[j-1]) / dx ### question 9 def resolution_lax(C , dt , dx , c_max , v_max): for i in range(1,len(C)): Q = debit(v_max , c_max ,C[i-1]) n = len(C[i]) for j in range(n): C_moy = (C[i - 1][j - 1] + C[i - 1][(j + 1) % n]) / 2 C[i][j] = C_moy - dt * (Q[(j + 1) % n] - Q[j - 1]) / (2 * dx) # mise en évidence des tracés def graphique(fonction, c1, c2, titre,ax1): C = [C_depart(La, dx, c1, c2, d1, d2)] t = dt while t < Temps: C.append([0]*len(C[0])) t += dt fonction(C, dt, dx, c_max, v_max) ax1.plot(np.arange(0,La,dx),C[0],linestyle='--') for i in range(1,len(C)-10,10): plt.plot(np.arange(0,La,dx),C[i]) ax1.plot(np.arange(0,La,dx),C[-1],linestyle=':') ax1.axvline(x=d1, linestyle='-.', color="grey") ax1.axvline(x=d2, linestyle='-.', color="grey") ax1.set_ylim([0, 0.2]) ax1.set_title(titre) La = 8500 # 8.5 km Temps= 100 # seconde dx= 50 # m dt= 1 # s d1 = La / 3 d2 = 2 * d1 fig = plt.figure() graphique(resolution, 0.01, 0.03, "Cas 1 Euler Avant", fig.add_subplot(321)) graphique(resolution, 0.12, 0.16, "Cas 2 Euler Avant", fig.add_subplot(322)) graphique(resolution_arriere, 0.01, 0.03, "Cas 1 Euler Arrière", fig.add_subplot(323)) graphique(resolution_arriere, 0.12, 0.16, "Cas 2 Euler Arrière", fig.add_subplot(324)) graphique(resolution_lax, 0.01, 0.03, "Cas 1 Lax-Friedriechs", fig.add_subplot(325)) graphique(resolution_lax, 0.12, 0.16, "Cas 2 Lax-Friedriechs", fig.add_subplot(326)) fig.tight_layout() plt.show()