Hello tout le monde !
Ecoutez les amis, je travaille actuellement sur un projet consistant à modéliser la trajectoire d'une balle de tennis...en prenant en compte les frottements et la rotation de la balle.
Du coup en projetant sur mes 3 axes j'obtient:
dvx/dt= -(h/m)*sqrt(vx**2+vy**2+vz**2)*vx+ (l/m)*(wy*vz-wz*vy)
dvy/dt= -(h/m)*sqrt(vx**2+vy**2+vz**2)*vy+ (l/m)*(wz*vx-wx*vz)-g
dvz/dt= -(h/m)*sqrt(vx**2+vy**2+vz**2)*vy+ (l/m)*(wx*vy-wy*vx)
Du coup j'ai programmé tout çasur python en m'aidant du module solve_ivp...sauf que je n'obtient rien….
Je demande donc de l'aide à un un bon samaritin pour trouver ce qui colle pas
Merci d'avance
Mon code:
Code:import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp from mpl_toolkits.mplot3d import Axes3D # global g, h, l, m, wx, wy,wz,x0,y0,z0,v0x,v0y,v0z g=9.81 #accélération de la pesanteur m=56.7E-03 #masse de la balle x0=0 #position du joueur sur le terrain y0=0 #hauteur initiale du lancé z0=2 v0=20 #vitesse du service θ=10 #angle de lancé par rapport à l'horizontal α=40 v0x=v0*np.cos(θ*np.pi/180) #vecteur vitesse projeté sur l'axe x v0z=v0*np.sin(θ*np.pi/180) #vecteur vitesse projeté sur l'axe x v0y=v0*np.sin(α*np.pi/180) p=1.225 #masse volumique de l'air R=3.3E-02 #rayon de la balle S=4*np.pi*R**2 C=0.4 #coeffcient de trainée h=0.5*p*C*S wx, wy, wz= 0,0,100 #vecteur rotation en rad/s l=0.17 def trajectoire(t, u): global g,l,m,wx,wy,wz,h x = u[0] y = u[1] z = u[2] vx = u[3] vy = u[4] vz = u[5] dudt = np.zeros(np.size(u)) # le système que l'on veut intégrer !' dudt[0] = u[3] # equation de la vitesse vx dudt[1] = u[4] # equation de la vitesse vy dudt[2] = u[5] dudt[2] = -(h/m)*np.sqrt(u[3]**2+u[4]**2+u[5]**2)*u[3]+(l/m)*(wy*u[5]-wz*u[4]) # equation de l'accélération sur x dudt[3] = -(h/m)*np.sqrt(u[3]**2+u[4]**2+u[5]**2)*u[4]+(l/m)*(wz*u[3]-wx*u[5]) # equation de l'accélération sur y dudt[4] = -(h/m)*np.sqrt(u[3]**2+u[4]**2+u[5]**2)*u[5]+(l/m)*(wx*u[4]-wy*u[3])-g return dudt tinit = 0 tfinal = 2.1 nt = 1001 t = np.linspace(tinit,tfinal,nt) # #************************************************************************ # method="RK45" "RK23" "Radau" "BDF" "LSODA" #************************************************************************ ## option # t_eval=tv , rtol = 1.0e-3 , atol = 1.0e-6 # ,rtol = 1.0e-6 , atol = 1.0e-8 # dense_output=True, # max_step = 0.01, events=None # résolution du système d'EDP sous forme d'EDO # conditions initiales Uinit = [0,0,2,27.63,4.8,27.36] # x0,y0, vx0, vy0 sol = solve_ivp(trajectoire, (tinit, tfinal), y0=Uinit, args=(wx,wy,wz,g,l,h,m) ) x = sol.y[0] y = sol.y[1] z = sol.y[2] vx = sol.y[3] vy = sol.y[4] vz = sol.y[5] i=0 for k in t: if z[i]<0: break i=i+1 fig, ax=plt.subplots() ax = plt.axes(projection='3d') ax.set_xlabel("Longueur (m)") ax.set_ylabel("Largeur (m)") ax.set_zlabel("Hauteur (m)") ax.plot(x[:i],y[:i],z[:i],lw=2) plt.show()
-----
