Skip to content
Extraits de code Groupes Projets
Valider 2c42fe0c rédigé par Lilian Vanderveken's avatar Lilian Vanderveken
Parcourir les fichiers

02/02/21

parent 7e4e1441
Branches master
Aucune étiquette associée trouvée
Aucune requête de fusion associée trouvée
Masquer les modifications d'espaces
En ligne Côte à côte
Affichage de
avec 0 ajout et 241 suppressions
TP1.ipynb supprimé 100644 → 0
+ 0
− 241
Voir le fichier @ 7e4e1441
%% Cell type:markdown id: tags:
# Oscillateur harmonique
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=6)
def HarmOsc(k,m,t,x0,mode):
dt=t[1]-t[0]
x=np.zeros((2,np.shape(t)[0]))
x[:,0]=x0
M=np.array([[0,-k],[1/m,0]])
if mode=='Euler avant':
for i in range(np.shape(x)[1]-1):
M=np.array([[0,-k],[1/m,0]])
x[:,i+1]=(np.identity(2)+dt*M).dot(x[:,i])
elif mode=='Euler arrière':
for i in range(np.shape(x)[1]-1):
M=np.array([[0,-k],[1/m,0]])
x[:,i+1]=np.linalg.solve((np.identity(2)-dt*M),x[:,i])
elif mode=='Euler centré':
for i in range(np.shape(x)[1]-1):
M=np.array([[0,-k],[1/m,0]])
x[:,i+1]=np.linalg.solve((np.identity(2)-dt/2*M),(np.identity(2)+dt/2*M).dot(x[:,i]))
elif mode=='Leapfrog':
x[:,1]=x[:,0]
for i in range (1,np.shape(x)[1]-1):
x[:,i+1]=x[:,i-1]+2*dt*M.dot(x[:,i])
elif mode=='Heun':
for i in range(np.shape(x)[1]-1):
xx=x[:,i]+dt*M.dot(x[:,i])
x[:,i+1]=x[:,i]+(dt/2)*(M.dot(x[:,i])+M.dot(xx))
elif mode=='RK4':
for i in range(np.shape(x)[1]-1):
k1=dt*M.dot(x[:,i])
k2=dt*M.dot(x[:,i]+k1/2)
k3=dt*M.dot(x[:,i]+k2/2)
k4=dt*M.dot(x[:,i]+k3)
x[:,i+1]=x[:,i]+(k1+2*k2+2*k3+k4)/6
elif mode=='Verlet':
for i in range(np.shape(x)[1]-1):
pp=x[0,i]-dt*k*x[1,i]
x[1,i+1]=x[1,i]+(dt/2)*(x[0,i]/m+pp/m)
x[0,i+1]=x[0,i]-(dt/2)*(x[1,i]+x[1,i+1])
elif mode=='Stormer_Verlet':
for i in range(np.shape(x)[1]-1):
pp=x[0,i]-(dt/2)*k*x[1,i]
x[1,i+1]=x[1,i]+(dt)*(pp/m)
x[0,i+1]=pp-(dt/2)*(x[1,i+1])
else:
print('Invalid integration method')
return(x)
def Energie_calc(x,k,m):
return((x[0,:]*x[0,:])/(2*m)+k*(x[1,:]*x[1,:])/2)
```
%% Cell type:code id: tags:
``` python
k=1.
m=1.0
omega=np.sqrt(k/m)
lmb=2*np.pi/omega
dt=0.001
T=50
numt=int(np.floor(T/dt))+1
t=np.linspace(0,T,numt)
#Conditions initiales
x0=[0,1]
xsol=np.zeros((2,numt))
xsol[1,:]=x0[1]*np.cos(omega*t)+x0[0]/(m*omega)*np.sin(omega*t)
xsol[0,:]=(-x0[1]*omega*np.sin(omega*t)+x0[0]/(m*omega)*omega*np.cos(omega*t))*m
x_EUav=HarmOsc(k,m,t,x0,'Euler avant')
x_EUar=HarmOsc(k,m,t,x0,'Euler arrière')
x_Verl=HarmOsc(k,m,t,x0,'Verlet')
x_StVer=HarmOsc(k,m,t,x0,'Stormer_Verlet')
x_Heun=HarmOsc(k,m,t,x0,'Heun')
x_RK4=HarmOsc(k,m,t,x0,'RK4')
x_EUcen=HarmOsc(k,m,t,x0,'Euler centré')
E_Verl=Energie_calc(x_Verl,k,m)
E_StVerl=Energie_calc(x_StVer,k,m)
E_RK4=Energie_calc(x_RK4,k,m)
E_Heun=Energie_calc(x_Heun,k,m)
E_EUcen=Energie_calc(x_EUcen,k,m)
E_EUar=Energie_calc(x_EUar,k,m)
E_EUav=Energie_calc(x_EUav,k,m)
Esol=Energie_calc(xsol,k,m)
```
%% Cell type:code id: tags:
``` python
%matplotlib widget
fig,ax=plt.subplots(figsize=(10,10))
plt.plot(t,E_Verl,label='Verlet')
plt.plot(t,E_StVerl,label='Stormer_Verlet')
plt.plot(t,E_Heun,label='Heun')
plt.plot(t,E_RK4,label='RK4')
plt.plot(t,E_EUcen,label='Euler centré',marker='D')
plt.plot(t,E_EUav,label='Euler avant')
plt.plot(t,E_EUar,label='Euler arrière')
plt.plot(t,Esol,label='Theoretical')
plt.ylim((0.45,0.525))
plt.xlim((-1,50))
plt.legend(loc='best',fontsize='small')
plt.xlabel('t')
plt.ylabel('Energie')
plt.title('Oscillateur harmonique: dt= %.5f'%(dt))
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
import pandas as pd
url="https://raw.githubusercontent.com/plotly/datasets/master/tips.csv"
df = pd.read_csv(url)
print(df)
df.plot(kind="line")
plt.scatter("total_bill", "tip",data=df)
plt.xlabel("Total Bill")
plt.ylabel("Tip")
plt.show()
```
%% Output
total_bill tip sex smoker day time size
0 16.99 1.01 Female No Sun Dinner 2
1 10.34 1.66 Male No Sun Dinner 3
2 21.01 3.50 Male No Sun Dinner 3
3 23.68 3.31 Male No Sun Dinner 2
4 24.59 3.61 Female No Sun Dinner 4
.. ... ... ... ... ... ... ...
239 29.03 5.92 Male No Sat Dinner 3
240 27.18 2.00 Female Yes Sat Dinner 2
241 22.67 2.00 Male Yes Sat Dinner 2
242 17.82 1.75 Male No Sat Dinner 2
243 18.78 3.00 Female No Thur Dinner 2
[244 rows x 7 columns]
%% Cell type:code id: tags:
``` python
!jupyter lab --version
```
%% Output
2.2.9
0% ou .
You are about to add 0 people to the discussion. Proceed with caution.
Terminez d'abord l'édition de ce message.
Veuillez vous inscrire ou vous pour commenter