TP1 Intro Modelisation

EA.py

@@ -0,0 +1,31 @@
+import numpy as np
+import sympy
+import matplotlib.pyplot as plt
+
+def EA1(xi,alpha):
+    return 1-np.sqrt(1-2*alpha*(1-alpha)*(1-np.cos(xi)))
+
+def EA2(xi,alpha):
+    return 1-np.sqrt(1-(alpha**2)*(1-alpha**2)*((1-np.cos(xi))**2))
+
+def EA3(xi,alpha):
+    theta=xi/2
+    return 1-np.sqrt((1-2*alpha*(np.sin(theta)**2)*(1-(1-alpha)*(np.cos(theta)**2)))**2
+        +(2*alpha*np.cos(theta)*np.sin(theta)*(1+(1-alpha)*(np.sin(theta)**2)))**2)
+
+xi = np.linspace(0,np.pi/8,1000)
+
+for alpha in [0.2,0.5,0.8]:
+    fig=plt.figure(1)
+    plt.title(r"Erreur d'amplitude en fonction de $\zeta$ pour CFL = "+str(alpha))
+    plt.xlabel(r'$\zeta$')
+    plt.ylabel(r'$E_A$')
+    plt.axis([-0.02 , np.pi/8+0.02 , 1e-15, 1e-1])
+    plt.semilogy(xi,EA1(xi,alpha),'-g',label="Schema 1")
+    plt.semilogy(xi,EA2(xi,alpha),'-b',label="Schema 2")
+    plt.semilogy(xi,EA3(xi,alpha),'-r',label="Schema 3")
+    plt.legend()
+    plt.grid(True)
+    plt.show(block=False)
+    #plt.savefig("EA"+"_CFL ="+str(alpha)+".png")
+    plt.close('all')

EP.py

@@ -0,0 +1,42 @@
+import numpy as np
+import math
+import matplotlib.pyplot as plt
+import EA
+
+def acos(X):
+    Y=np.zeros_like(X)
+    for i in range(len(X)):
+        Y[i]=math.acos(X[i])
+    return Y
+
+def EP1(xi,alpha):
+    return alpha*xi - acos((1-alpha*(1-np.cos(xi)))
+                           /(1-EA.EA1(xi,alpha)))
+
+def EP2(xi,alpha):
+    return alpha*xi - acos((1-(alpha**2)*(1-np.cos(xi)))
+                            /(1-EA.EA2(xi,alpha)))
+
+def EP3(xi,alpha):
+    theta=xi/2
+    return alpha*xi - acos((1-2*alpha*(np.sin(theta)**2)*(1-(1-alpha)*(np.cos(theta)**2)))
+                            /(1-EA.EA3(xi,alpha)))
+
+xi = np.linspace(0,np.pi/8,1000)
+
+#print(EP1(xi,0.5),"\n",EP3(xi,0.5))
+
+for alpha in [0.2,0.5,0.8]:
+    fig=plt.figure(1)
+    plt.title(r"Erreur de phase en fonction de $\zeta$ pour CFL = "+str(alpha))
+    plt.xlabel(r'$\zeta$')
+    plt.ylabel(r'$E_\phi$')
+    plt.axis([-0.02 , np.pi/8+0.02 , -0.0010, 0.0030])
+    plt.plot(xi,abs(EP1(xi,alpha)),'-g',label="Schema 1")
+    plt.plot(xi,abs(EP2(xi,alpha)),'-b',label="Schema 2")
+    plt.plot(xi,abs(EP3(xi,alpha)),'-r',label="Schema 3")
+    plt.legend(loc='upper left')
+    plt.grid(True)
+    plt.show(block=False)
+    #plt.savefig("EP"+"_CFL ="+str(alpha)+".png")
+    plt.close('all')

tp1.py

@@ -0,0 +1,52 @@
+# -*- coding: cp1252 -*-
+import numpy as np
+import matplotlib.pyplot as plt
+import tp1_lib as tp
+
+## Initialisation des données 
+Lx = 1.0                     # Demi longueur de l'intervalle       
+a = 2.0                      # Vitesse du signal 
+tmax = 2.0                   # Temps max de la simulation (=T)
+
+J= 201                       # Nombre de noeuds
+X = np.linspace(-Lx,Lx,J)    # Abscisses des noeuds
+dx = X[1]-X[0]               # Pas de maillage
+cfl = 0.8
+#cfl = 0.5
+#cfl = 0.2                    # Nombre CFL 
+dt = cfl*dx/a                # Valeur du pas de temps 
+
+## Initialisation de la solution 
+Uini = tp.Uinit(X,J)         
+U1 = Uini   # Schema 1
+U2 = Uini   # Schema 2
+U3 = Uini   # Schema 3
+
+for tmax in (1,5,10):
+    ##Boucle en temps (tant que time < tmax) 
+    time = 0.0                     
+    while time < tmax :
+        dt_reel = min(dt, tmax-time)   # dt_reel = dt sauf si tmax - time < dt 
+        alpha = a*dt_reel/dx           # Nombre CFL effectif 
+        print ('TIME ITERATION : ' , time , ' \ ' , tmax)
+        U1 = tp.schema1(U1,J,alpha)    # Mise à jour de la solution schema 1
+        U2 = tp.schema2(U2,J,alpha)    # Mise à jour de la solution schema 2
+        U3 = tp.schema3(U3,J,alpha)    # Mise à jour de la solution schema 3
+        time += dt_reel                # Incrémente le temps 
+     
+    ## Affichage des résultats
+    fig=plt.figure(1)
+    plt.title("Solution u en fonction de x pour CFL = "+str(cfl)+" et T = "+str(tmax))
+    plt.xlabel('x')
+    plt.ylabel('u')
+    plt.axis([-1 , 1 , -0.25, 1.25])
+    plt.plot(X,U1,'-g',X,U2,'-b',X,U3,'-r',X,Uini,'-k')
+    plt.legend(['Schema1','Schema2', 'Schema3', 'Exacte'],loc='best')
+    plt.grid(True)     
+    plt.show(block=False)
+    plt.pause(10.0)  ## Pause pour voir la figure avant sauvegarde
+    #plt.savefig("u=X(x)"+"T="+str(tmax)+"_CFL ="+str(cfl)+".png")  ##sauve figure
+    plt.close('all')
+
+
+

tp1_lib.py

@@ -0,0 +1,34 @@
+import numpy as np
+
+def khi(X):
+    khi = np.zeros_like(X)
+    khi[(X >= -1/4) & (X <= 1/4)] = 1
+    return khi
+
+def Uinit(X,J):
+    Uinit = np.zeros(J)
+    Uinit[:] = 256 * ((X[:]-1/4)**2) * ((X[:]+1/4)**2) * khi(X)[:]
+    return Uinit
+
+"""def Uinit(X,J):
+    Uinit = np.zeros(J)
+    Uinit[:] = khi(X)[:]
+    return Uinit"""
+
+def schema1(U0,J,alpha):
+    U0=np.append(U0[J-2],U0)
+    U1=U0[1:J+1] - alpha*(U0[1:J+1]-U0[0:J])
+    return U1
+
+def schema2(U0,J,alpha):
+    U0=np.append(U0[J-2],U0)
+    U0=np.append(U0,U0[1])
+    U1=U0[1:J+1]-(alpha/2)*(U0[2:J+2]-U0[0:J])+(0.5*alpha**2)*(U0[2:J+2]-2*U0[1:J+1]+U0[0:J])
+    return U1
+
+def schema3(U0,J,alpha):
+    U0=np.append(U0[J-2],U0)
+    U0=np.append(U0[J-3],U0)
+    U0=np.append(U0,U0[1])
+    U1=U0[2:J+2]-alpha*(U0[2:J+2]-U0[1:J+1])-((alpha/4)*(1-alpha)*(U0[3:J+3]-U0[2:J+2]-U0[1:J+1]+U0[0:J]))             
+    return U1