Simulación de una pelota que rebota

$$m\frac{d\vec r}{dt} = -m g \hat j$$

tomemos en cuenta que $v_x = \dot x$ y $v_y = \dot y$ y con ello obtenemos: $$\begin{array}\\ \dot x = v_x\\ \dot v_x = 0\\ \dot y = v_y\\ \dot v_y = -g \end{array}$$

from pylab import *
from scipy.integrate import *
from matplotlib import animation,rc 
from IPython.display import HTML
%config InlineBackend.figure_format = 'retina' # Solo para pantallas de alta resolucion

def tiro_sf(x_v,t,g=9.8):
    # x_v son       x,v_x,y,v_y
    # con indices   0, 1 ,2, 3
    dx = x_v[1]
    dvx = 0
    dy = x_v[3]
    dvy = -g
    return array([dx,dvx,dy,dvy])

t = linspace(0,0.5)
cond_ini = [0,2,0,2]
x_v = odeint(tiro_sf,cond_ini,t)
plot(x_v[:,0],x_v[:,2],'o-')
grid()
show()

def trayectoria(sist_ecs,ci,dt=0.01,N=50):
    x_t,y_t = zeros(N),zeros(N)
    x_t[0],y_t[0] = ci[0],ci[2]
    for i in range(N-1):
        x_v = odeint(sist_ecs,ci,[0,dt])
        ci = x_v[-1,:]
        if x_v[0,2]*x_v[1,2] < 0:
            print("Hay un rebote, con indice: ",i)
            ci = x_v[0,:]
            ci = rebote(sist_ecs,ci,dt)
            ci[3] = -ci[3]
        x_t[i+1],y_t[i+1] = ci[0],ci[2]
    return x_t,y_t

def rebote(sist_ecs,ci,dt,error=1e-6,max_iter=100):
    # calcula el rebote unicamente en y, es decir con el indice 2
    dif,i = 1,0
    while dif > error or i < max_iter:
        dt = 0.5*dt
        x_v_rebotada = odeint(sist_ecs,ci,[0,dt])
        if x_v_rebotada[0,2]*x_v_rebotada[1,2] < 0:
            ci = x_v_rebotada[0,:]
        else:
            ci = x_v_rebotada[1,:]
        dif = abs( ci[2] )
        i += 1
    return ci

x,y = trayectoria(tiro_sf,[0,2,0,2],N=100)
plot(x,y,'o-')
grid()
show()

Hay un rebote, con indice:  40
Hay un rebote, con indice:  81

print(y[40],y[41],y[42])

0.01599999941775676 9.087154983575799e-33 0.0195099999566561

$$m\frac{d\vec r}{dt} = -m g \hat j + m\gamma v \hat v$$

tomemos en cuenta que $v_x = \dot x$ y $v_y = \dot y$ y con ello obtenemos: $$\begin{array}\\ \dot x = v_x\\ \dot v_x = -\gamma v_x\\ \dot y = v_y\\ \dot v_y = -g -\gamma v_y \end{array}$$

def tiro_cf_pv(x_v,t,g=9.8,gam=0.1):
    # x_v son       x,v_x,y,v_y
    # con indices   0, 1 ,2, 3
    dx = x_v[1]
    dvx = - gam*x_v[1]
    dy = x_v[3]
    dvy = -g -gam*x_v[3]
    return array([dx,dvx,dy,dvy])

x_fpv,y_fpv = trayectoria(tiro_cf_pv,[0,2,0,2],N=100)
plot(x,y,'o-',label="Sin fricción")
plot(x_fpv,y_fpv,'o-',label="Fricción $\sim v$")
legend(),grid()
show()

Hay un rebote, con indice:  40
Hay un rebote, con indice:  81

$$m\frac{d\vec r}{dt} = -m g \hat j + m\gamma v^2 \hat v$$

tomemos en cuenta que $v_x = \dot x$ y $v_y = \dot y$ y con ello obtenemos: $$\begin{array}\\ \dot x = v_x\\ \dot v_x = -\gamma v v_x\\ \dot y = v_y\\ \dot v_y = -g -\gamma v v_y \end{array}$$ donde $v=\sqrt{v_x ^2 + v_y ^2}$

def tiro_cf_pv2(x_v,t,g=9.8,gam=0.1):
    # x_v son       x,v_x,y,v_y
    # con indices   0, 1 ,2, 3
    v = sqrt(x_v[1]**2 + x_v[3]**2)
    dx = x_v[1]
    dvx = - gam*v*x_v[1]
    dy = x_v[3]
    dvy = -g -gam*v*x_v[3]
    return array([dx,dvx,dy,dvy])

x_fpv2,y_fpv2 = trayectoria(tiro_cf_pv2,[0,2,0,2],N=100)
plot(x,y,'-',label="Sin fricción")
plot(x_fpv,y_fpv,'-',label="Fricción $\sim v$")
plot(x_fpv2,y_fpv2,'-',label="Fricción $\sim v^2$")
legend(),grid()
show()

Hay un rebote, con indice:  40
Hay un rebote, con indice:  79

$$m\frac{d\vec r}{dt} = -m g \hat j + m\gamma v^2 \hat v - mf_0\hat i$$

tomemos en cuenta que $v_x = \dot x$ y $v_y = \dot y$ y con ello obtenemos: $$\begin{array}\\ \dot x = v_x\\ \dot v_x = -\gamma v v_x - f_0\\ \dot y = v_y\\ \dot v_y = -g -\gamma v v_y \end{array}$$ donde $v=\sqrt{v_x ^2 + v_y ^2}$

def tiro_cv(x_v,t,g=9.8,gam=0.1,f0=3):
    # x_v son       x,v_x,y,v_y
    # con indices   0, 1 ,2, 3
    v = sqrt(x_v[1]**2 + x_v[3]**2)
    dx = x_v[1]
    dvx = - gam*v*x_v[1] - f0
    dy = x_v[3]
    dvy = -g -gam*v*x_v[3]
    return array([dx,dvx,dy,dvy])

x_cv,y_cv = trayectoria(tiro_cv,[0,2,0,2],N=100)
plot(x,y,'-',label="Sin fricción")
plot(x_fpv,y_fpv,'-',label="Fricción $\sim v$")
plot(x_fpv2,y_fpv2,'-',label="Fricción $\sim v^2$")
plot(x_cv,y_cv,'-',label="Fricción $\sim v^2$ y $f_0$")
legend(),grid()
show()

Hay un rebote, con indice:  40
Hay un rebote, con indice:  80

%%capture

fig,ejes = plt.subplots(figsize=(3.5,2),dpi=200)
ejes.set_xlim((0,2)),ejes.set_ylim((0,.21)),ejes.grid()
ejes.set_title(r'Tiro parabólico')

t1, = ejes.plot([],[],'-',color='blue')
t2, = ejes.plot([],[],'-',color='orange')
t3, = ejes.plot([],[],'-',color='green')
t4, = ejes.plot([],[],'-',color='red')

def animar(i):
    t1.set_data([x[:i]],[y[:i]])
    t2.set_data([x_fpv[:i]],[y_fpv[:i]])
    t3.set_data([x_fpv2[:i]],[y_fpv2[:i]])
    t4.set_data([x_cv[:i]],[y_cv[:i]])
    return (t1,t2,t3,t4,)

pelicula = animation.FuncAnimation(fig,animar,frames=100,interval=40,blit=True)

HTML( pelicula.to_html5_video() )