commit 5473ade7817b86bfad632acc90b821db9d2eb5d6 Author: mbui Date: Tue May 9 22:19:36 2023 +0200 TP1 Intro Modelisation diff --git a/EA.py b/EA.py new file mode 100644 index 0000000..68cd8cb --- /dev/null +++ b/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') diff --git a/EP.py b/EP.py new file mode 100644 index 0000000..da90ab3 --- /dev/null +++ b/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') \ No newline at end of file diff --git a/Fig/EA_CFL =0.2.png b/Fig/EA_CFL =0.2.png new file mode 100644 index 0000000..cd66bd5 Binary files /dev/null and b/Fig/EA_CFL =0.2.png differ diff --git a/Fig/EA_CFL =0.5.png b/Fig/EA_CFL =0.5.png new file mode 100644 index 0000000..4f58cce Binary files /dev/null and b/Fig/EA_CFL =0.5.png differ diff --git a/Fig/EA_CFL =0.8.png b/Fig/EA_CFL =0.8.png new file mode 100644 index 0000000..1a07a36 Binary files /dev/null and b/Fig/EA_CFL =0.8.png differ diff --git a/Fig/EP_CFL =0.2.png b/Fig/EP_CFL =0.2.png new file mode 100644 index 0000000..c4f9a35 Binary files /dev/null and b/Fig/EP_CFL =0.2.png differ diff --git a/Fig/EP_CFL =0.5.png b/Fig/EP_CFL =0.5.png new file mode 100644 index 0000000..8d9de9f Binary files /dev/null and b/Fig/EP_CFL =0.5.png differ diff --git a/Fig/EP_CFL =0.8.png b/Fig/EP_CFL =0.8.png new file mode 100644 index 0000000..d24487b Binary files /dev/null and b/Fig/EP_CFL =0.8.png differ diff --git a/Fig/u=X(x)T=10_CFL =0.8.png b/Fig/u=X(x)T=10_CFL =0.8.png new file mode 100644 index 0000000..8553bec Binary files /dev/null and b/Fig/u=X(x)T=10_CFL =0.8.png differ diff --git a/Fig/u=X(x)T=1_CFL =0.8.png b/Fig/u=X(x)T=1_CFL =0.8.png new file mode 100644 index 0000000..384def4 Binary files /dev/null and b/Fig/u=X(x)T=1_CFL =0.8.png differ diff --git a/Fig/u=X(x)T=5_CFL =0.8.png b/Fig/u=X(x)T=5_CFL =0.8.png new file mode 100644 index 0000000..fb1a125 Binary files /dev/null and b/Fig/u=X(x)T=5_CFL =0.8.png differ diff --git a/Fig/u=f(x)CFL =0.2.png b/Fig/u=f(x)CFL =0.2.png new file mode 100644 index 0000000..fac55ce Binary files /dev/null and b/Fig/u=f(x)CFL =0.2.png differ diff --git a/Fig/u=f(x)CFL =0.5.png b/Fig/u=f(x)CFL =0.5.png new file mode 100644 index 0000000..b4453a9 Binary files /dev/null and b/Fig/u=f(x)CFL =0.5.png differ diff --git a/Fig/u=f(x)CFL =0.8.png b/Fig/u=f(x)CFL =0.8.png new file mode 100644 index 0000000..3f48a41 Binary files /dev/null and b/Fig/u=f(x)CFL =0.8.png differ diff --git a/Fig/u=f(x)T=10_CFL =0.2.png b/Fig/u=f(x)T=10_CFL =0.2.png new file mode 100644 index 0000000..ed320c8 Binary files /dev/null and b/Fig/u=f(x)T=10_CFL =0.2.png differ diff --git a/Fig/u=f(x)T=10_CFL =0.5.png b/Fig/u=f(x)T=10_CFL =0.5.png new file mode 100644 index 0000000..22a3fcb Binary files /dev/null and b/Fig/u=f(x)T=10_CFL =0.5.png differ diff --git a/Fig/u=f(x)T=10_CFL =0.8.png b/Fig/u=f(x)T=10_CFL =0.8.png new file mode 100644 index 0000000..4796c96 Binary files /dev/null and b/Fig/u=f(x)T=10_CFL =0.8.png differ diff --git a/Fig/u=f(x)T=1_CFL =0.2.png b/Fig/u=f(x)T=1_CFL =0.2.png new file mode 100644 index 0000000..b2d8600 Binary files /dev/null and b/Fig/u=f(x)T=1_CFL =0.2.png differ diff --git a/Fig/u=f(x)T=1_CFL =0.5.png b/Fig/u=f(x)T=1_CFL =0.5.png new file mode 100644 index 0000000..3d44042 Binary files /dev/null and b/Fig/u=f(x)T=1_CFL =0.5.png differ diff --git a/Fig/u=f(x)T=1_CFL =0.8.png b/Fig/u=f(x)T=1_CFL =0.8.png new file mode 100644 index 0000000..6c05338 Binary files /dev/null and b/Fig/u=f(x)T=1_CFL =0.8.png differ diff --git a/Fig/u=f(x)T=5_CFL =0.2.png b/Fig/u=f(x)T=5_CFL =0.2.png new file mode 100644 index 0000000..08ec326 Binary files /dev/null and b/Fig/u=f(x)T=5_CFL =0.2.png differ diff --git a/Fig/u=f(x)T=5_CFL =0.5.png b/Fig/u=f(x)T=5_CFL =0.5.png new file mode 100644 index 0000000..9b3aedd Binary files /dev/null and b/Fig/u=f(x)T=5_CFL =0.5.png differ diff --git a/Fig/u=f(x)T=5_CFL =0.8.png b/Fig/u=f(x)T=5_CFL =0.8.png new file mode 100644 index 0000000..55d2867 Binary files /dev/null and b/Fig/u=f(x)T=5_CFL =0.8.png differ diff --git a/tp1.py b/tp1.py new file mode 100644 index 0000000..d2a699e --- /dev/null +++ b/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') + + + diff --git a/tp1_lib.py b/tp1_lib.py new file mode 100644 index 0000000..d1feac8 --- /dev/null +++ b/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