¿Cómo vibra un tambor cuando lo golpeas?

Analizar el problema de la membrana vibrante permite entender el funcionamiento de instrumentos de percusión tales como los tambores, timbales e incluso sistemas biológicos como el tímpano.

Referencias:

Considere un tambor (membrana) de radio $a$, entonces la función de onda en $\mathbb{R}^2$ para este sistema se puede escribir como, $$ \frac{1}{v}\frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} $$

donde $u\equiv u(x,y,t)$ es el desplazamiento transversal(elevación) y $v$ es la rapidez de propagación de la onda.

La forma habitual de encontrar la solución a esta ecuación es primero hacer un cambio de coordenadas, de cartesianas a polares y posteriormente considerar el método de separación de variable, tal que

$$ u(r, \theta, t) = R(r) \Theta(\theta) T (t)$$

Esta sustitución da como resultado tres ecuaciones diferenciales, una para cada variable de separación. Y cuyas soluciones nos ayudan a escribir los modos normales.

$$u_{nk}(r,\theta, t) = J_{n}(\lambda_{nk} r)(a_{nk}\cos{n\theta} + b_{nk}\sin{n\theta})\cos{(v\lambda_{nk}t)}$$$$u^{*}_{nk}(r,\theta, t) = J_{n}(\lambda_{nk}r)(a^{*}_{nk}\cos{n\theta} + b^{*}_{nk}\sin{n\theta})\sin{(v\lambda_{nk}t)})$$

para $n = 0,1,2,\dots$, $k = 1,2,3,\dots$, donde $J_{n}$ es la función de Bessel de orden $n$ de primera clase. Además, $$\lambda_{nk} = \frac{\alpha_{nk}}{a}$$

donde $\alpha_{nk}$ es el k-ésimo cero de $J_{n}(\lambda a)=0$. Esto es consecuencia de que $u$ sea cero en la frontera de la membrana, $r = a$.

Los coeficientes $a_{nk} , b_{nk}, a^{*}_{nk}$ y $b^{*}_{nk}$ se determinan de tal forma que cumplan con las condiciones iniciales:

$$u(r,\theta, 0) = f(r,\theta)$$$$u_{t}(r,\theta, 0) = g(r,\theta)$$

donde el primer termino es la geometría inicial y el segundo la rapidez inicial. Se puede demostrar que las expresiones para estos coeficientes se pueden escribir como:

\begin{align} a_{0k} &= \frac{1}{\pi a^2 J_{1}^{2}(\alpha_{0k})}\int_{0}^{2\pi}\int_{0}^{a}\; f(r,\theta)\, J_{0}(\lambda_{0k}r)\, r \, dr \, d\theta\\ a_{nk} &= \frac{2}{\pi a^2 J_{n+1}^{2}(\alpha_{nk})}\int_{0}^{2\pi}\int_{0}^{a}\; f(r,\theta)\, J_{n}(\lambda_{nk}r)\cos(n\theta)\, r \, dr \, d\theta\\ b_{nk} &= \frac{2}{\pi a^2 J_{n+1}^{2}(\alpha_{nk})}\int_{0}^{2\pi}\int_{0}^{a}\; f(r,\theta)\, J_{n}(\lambda_{nk}r)\sin(n\theta)\, r \, dr \, d\theta \end{align}

Y similarmente,

\begin{align} a^{*}_{0k} &= \frac{1}{\pi \,v\, \alpha_{0k}\,a J_{1}^{2}(\alpha_{0k})}\int_{0}^{2\pi}\int_{0}^{a}\; g(r,\theta)\, J_{0}(\lambda_{0k}r)\, r \, dr \, d\theta\\ a^{*}_{nk} &= \frac{2}{\pi\, v\,\alpha_{0k}\, a J_{n+1}^{2}(\alpha_{nk})}\int_{0}^{2\pi}\int_{0}^{a}\; g(r,\theta)\, J_{n}(\lambda_{nk}r)\cos(n\theta)\, r \, dr \, d\theta\\ b^{*}_{nk} &= \frac{2}{\pi\, v\,\alpha_{0k}\, a J_{n+1}^{2}(\alpha_{nk})}\int_{0}^{2\pi}\int_{0}^{a}\; g(r,\theta)\, J_{n}(\lambda_{nk}r)\sin(n\theta)\, r \, dr \, d\theta \end{align}

Estamos familiarizados con la función coseno, pero no tanto con la función de Bessel. Entonces, nuestra primera actividad será conocer su comportamiento.

\begin{align} u(r,\theta, t) &= \sum_{n=0}^{\infty}\sum_{k = 1}^{\infty}J_{n}(\lambda_{nk} r)(a_{nk}\cos{n\theta} + b_{nk}\sin{n\theta})\cos{(v\lambda_{nk}t)}\\ &+ \sum_{n=0}^{\infty}\sum_{k = 1}^{\infty}J_{n}(\lambda_{nk}r)(a^{*}_{nk}\cos{n\theta} + b^{*}_{nk}\sin{n\theta})\sin{(v\lambda_{nk}t)}) \end{align}



In [15]:

    
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import special
import numpy as np
from ipywidgets import *



In [16]:

    
r = np.linspace(0, 10,100)
for n in range(5):
    plt.plot(r, special.jn(n, r), label = '$J_{%s}(r)$'%n)
plt.xlabel('$r$', fontsize = 18)
plt.ylabel('$J_{n}(r)$', fontsize = 18)
plt.axhline(y = 0, color = 'k')
plt.legend(loc='center left', bbox_to_anchor=(1.0, 0.5), prop={'size': 14})
plt.show()

Por simplicidad vamos a suponer que $r = 1$ y determinar los ceros, significa encontrar todas las intersecciones de las curvas anteriores con el eje horizontal.

Ejemplo: Caso radialmente simétrico

Suponga que $a = 1$, $v = 1$ y que las condiciones iniciales son:

$$ f(r,\theta) = 1- r^4\quad\quad g(r,\theta) = 0$$

Dado que la rapidez inicial es cero, entonces $a^{*}_{nk} = b^{*}_{nk} = 0$



In [263]:

    
def f_shape(r):
    return 1 - r**4



In [264]:

    
r = np.r_[0:1:100j]
angle = np.r_[0:2*np.pi:200j]
r_shape = f_shape(r)



In [265]:

    
u = np.array([np.full(len(angle), radi) for radi in r_shape])
x = np.array([var_r * np.cos(angle) for var_r in r])
y = np.array([var_r * np.sin(angle) for var_r in r])



In [266]:

    
plt.figure(figsize = (6, 5))
plt.pcolor(x, y, u_good, cmap = 'viridis')
plt.axis('off')
plt.colorbar()
plt.show()

Y la solución para el desplazamiento en el tiempo es simplemente, \begin{equation} u(r,\theta, t) = \sum_{n=0}^{\infty}\sum_{k = 1}^{\infty}J_{n}(\lambda_{nk} r)(a_{nk}\cos{n\theta} + b_{nk}\sin{n\theta})\cos{(v\lambda_{nk}t)} \end{equation}

Entonces, solo será necesario encontrar $a_{nk}$ y $b_{nk}$.

Para resolver estas integrales haremos uso de sympy. Iniciemos con $a_{nk}$.

Primer caso $n>0$



In [130]:

    
from sympy import init_printing; init_printing(use_latex='mathjax')
import sympy as sym



In [131]:

    
r, theta, k  = sym.symbols('r theta k')
r, theta, k









    Out[131]:





$$\left ( r, \quad \theta, \quad k\right )$$



In [132]:

    
n = sym.Symbol('n', positive = True, integer=True)
n









    Out[132]:





$$n$$



In [133]:

    
def lamb(n,k):
    return sym.Symbol('lambda_%s%s'%(n,k), positive = True)



In [134]:

    
lamb(0,k)









    Out[134]:





$$\lambda_{0k}$$



In [135]:

    
f = 1 - r**4; f









    Out[135]:





$$- r^{4} + 1$$



In [136]:

    
integrand = f * sym.besselj(n, lamb(n,k) * r) * sym.cos(n *theta) * r
integrand









    Out[136]:





$$r \left(- r^{4} + 1\right) \cos{\left (n \theta \right )} J_{n}\left(\lambda_{nk} r\right)$$



In [137]:

    
ank = sym.Integral(integrand, (r, 0, 1), (theta, 0, 2*sym.pi))
ank









    Out[137]:





$$\int_{0}^{2 \pi}\int_{0}^{1} r \left(- r^{4} + 1\right) \cos{\left (n \theta \right )} J_{n}\left(\lambda_{nk} r\right)\, dr\, d\theta$$



In [138]:

    
solution = ank.doit()
solution









    Out[138]:





$$0$$

Entonces para cualquier $n>0$ no se tiene contribución.

Segundo caso $n=0$



In [139]:

    
integ = lambda n: f * sym.besselj(n, lamb(n,k) * r) * sym.cos(n*theta) * r



In [140]:

    
integ(0)









    Out[140]:





$$r \left(- r^{4} + 1\right) J_{0}\left(\lambda_{0k} r\right)$$



In [141]:

    
a0k = sym.Integral(integ(0), (r, 0, 1), (theta, 0, 2*sym.pi))
a0k









    Out[141]:





$$\int_{0}^{2 \pi}\int_{0}^{1} r \left(- r^{4} + 1\right) J_{0}\left(\lambda_{0k} r\right)\, dr\, d\theta$$



In [142]:

    
a0k_solution = a0k.doit()
a0k_solution









    Out[142]:





$$2 \pi \left(- \frac{4 J_{0}\left(\lambda_{0k}\right)}{\lambda_{0k}^{2}} + \frac{16 J_{1}\left(\lambda_{0k}\right)}{\lambda_{0k}^{3}} + \frac{32 J_{0}\left(\lambda_{0k}\right)}{\lambda_{0k}^{4}} - \frac{64 J_{1}\left(\lambda_{0k}\right)}{\lambda_{0k}^{5}}\right)$$



In [143]:

    
a0k_sol = 1/(sym.pi*sym.besselj(1, lamb(0,k))**2)*a0k_solution
a0k_sol









    Out[143]:





$$\frac{1}{J^{2}_{1}\left(\lambda_{0k}\right)} \left(- \frac{8 J_{0}\left(\lambda_{0k}\right)}{\lambda_{0k}^{2}} + \frac{32 J_{1}\left(\lambda_{0k}\right)}{\lambda_{0k}^{3}} + \frac{64 J_{0}\left(\lambda_{0k}\right)}{\lambda_{0k}^{4}} - \frac{128 J_{1}\left(\lambda_{0k}\right)}{\lambda_{0k}^{5}}\right)$$



In [144]:

    
sym.simplify(a0k_sol)









    Out[144]:





$$\frac{1}{\lambda_{0k}^{5} J^{2}_{1}\left(\lambda_{0k}\right)} \left(- 8 \lambda_{0k}^{3} J_{0}\left(\lambda_{0k}\right) + 32 \lambda_{0k}^{2} J_{1}\left(\lambda_{0k}\right) + 64 \lambda_{0k} J_{0}\left(\lambda_{0k}\right) - 128 J_{1}\left(\lambda_{0k}\right)\right)$$

Y para $b_{nk}$



In [109]:

    
integrand_b = f * sym.besselj(n, lamb(n,k) * r) * sym.sin(n *theta) * r
integrand_b









    Out[109]:





$$r \left(- r^{2} + 1\right) \sin{\left (n \theta \right )} J_{n}\left(\lambda_{nk} r\right)$$



In [111]:

    
bnk = sym.Integral(integrand_b, (r, 0, 1), (theta, 0, 2*sym.pi))
bnk









    Out[111]:





$$\int_{0}^{2 \pi}\int_{0}^{1} r \left(- r^{2} + 1\right) \sin{\left (n \theta \right )} J_{n}\left(\lambda_{nk} r\right)\, dr\, d\theta$$



In [112]:

    
solution_b = bnk.doit()
solution_b









    Out[112]:





$$0$$

¿Qué sucede para $n = 0$?



In [113]:

    
integ_b = lambda n: f * sym.besselj(n, lamb(n,k) * r) * sym.sin(n*theta) * r



In [114]:

    
integ_b(0)









    Out[114]:





$$0$$

Ahhh! el integrando es cero, entonces este termino también es cero.

\begin{equation} u(r,\theta, t) = \sum_{k = 1}^{\infty} a_{0k}J_{0}(\lambda_{0k} r)\cos{(v\lambda_{0k}t)} \end{equation}



In [145]:

    
a0k_sol









    Out[145]:





$$\frac{1}{J^{2}_{1}\left(\lambda_{0k}\right)} \left(- \frac{8 J_{0}\left(\lambda_{0k}\right)}{\lambda_{0k}^{2}} + \frac{32 J_{1}\left(\lambda_{0k}\right)}{\lambda_{0k}^{3}} + \frac{64 J_{0}\left(\lambda_{0k}\right)}{\lambda_{0k}^{4}} - \frac{128 J_{1}\left(\lambda_{0k}\right)}{\lambda_{0k}^{5}}\right)$$

Primero vamos a programar para algún modo $k$.



In [161]:

    
def a0k_sym(lambd):
    solucion = (-8*special.jn(0, lambd)/lambd**2 
    +32*special.jn(1, lambd)/lambd**3 + 
     64*special.jn(0, lambd)/lambd**4 -
        128*special.jn(1, lambd)/lambd**5)/special.jn(1, lambd)**2
    return solucion



In [208]:

    
def tambor(v, kth_zero, nt, t):
    r = np.r_[0:1:100j]
    angle = np.r_[0:2*np.pi:200j]
    ceros = special.jn_zeros(0, nt)
    lambd = ceros[kth_zero]
    u_r = a0k_sym(lambd)*special.jn(0, lambd * r) * np.cos(lambd * v * t)
    u = np.array([np.full(len(angle), u_rs) for u_rs in u_r])
    x = np.array([var_r * np.cos(angle) for var_r in r])
    y = np.array([var_r * np.sin(angle) for var_r in r])
    return x, y, u



In [243]:

    
x1, y1, u1 = tambor(1, 0, 15, 7)



In [244]:

    
plt.figure(figsize = (6, 5))
plt.pcolor(x1 , y1 , u1, cmap = 'viridis')
plt.axis('off')
plt.colorbar()
plt.show()



In [246]:

    
def tambor_nk(t = 0, kth=0):
    fig = plt.figure(figsize = (6,5))
    ax = fig.add_subplot(1, 1, 1)
    x, y, u = tambor(1, kth, 50, t)
    im = ax.pcolor(x, y, u, cmap = 'viridis')
    ax.set_xlim(xmin = -1.1, xmax = 1.1)
    ax.set_ylim(ymin = -1.1, ymax = 1.1)
    plt.colorbar(im)
    plt.axis('off')
    fig.canvas.draw()
interact_manual(tambor_nk, t = (0, 15,.01), n = (0, 10, 1), kth = (0, 10, 1));



In [249]:

    
def tambor_nk(t = 0, kth=0):
    fig = plt.figure(figsize = (6,5))
    ax = fig.add_subplot(1, 1, 1)
    x, y, u = tambor(1, kth, 50, t)
    im = ax.pcolor(x, y, u, cmap = 'viridis', vmax = 1.2, vmin = -1.2)
    ax.set_xlim(xmin = -1.1, xmax = 1.1)
    ax.set_ylim(ymin = -1.1, ymax = 1.1)
    plt.colorbar(im)
    plt.axis('off')
    fig.canvas.draw()
interact(tambor_nk, t = (0, 15,.01), n = (0, 10, 1), kth = (0, 10, 1));

Y ahoara, la solución completa.



In [256]:

    
def tambor_n_allk(v, nk_zeros, t):
    r = np.r_[0:1:100j]
    angle = np.r_[0:2*np.pi:200j]
    ceros = special.jn_zeros(0, nk_zeros)
    lambd = ceros[0]
    u_r = a0k_sym(lambd)*special.jn(0, lambd * r) * np.cos(lambd * v * t)
    u0 = np.array([np.full(len(angle), u_rs) for u_rs in u_r])
    for cero in range(1, nk_zeros):
        lambd = ceros[cero]
        u_r = a0k_sym(lambd)*special.jn(0, lambd * r) * np.cos(lambd * v * t)
        u = np.array([np.full(len(angle), u_rs) for u_rs in u_r])
        u0 += u
    x = np.array([var_r * np.cos(angle) for var_r in r])
    y = np.array([var_r * np.sin(angle) for var_r in r])
    return x, y, u0



In [257]:

    
def tambor_0(t = 0):
    fig = plt.figure(figsize = (6,5))
    ax = fig.add_subplot(1, 1, 1)
    x, y, u = tambor_n_allk(1, 15, t)
    im = ax.pcolor(x, y, u, cmap = 'viridis')
    ax.set_xlim(xmin = -1.1, xmax = 1.1)
    ax.set_ylim(ymin = -1.1, ymax = 1.1)
    plt.colorbar(im)
    plt.axis('off')
    fig.canvas.draw()
interact_manual(tambor_0, t = (0, 15,.01));

Fíjise bien, la condición inicial en $t = 0$, se cumple para la solución encontrada.

Tarea

Problema 1. Suponga que $a = 1$, $v = 1$ y que las condiciones iniciales son:

$$ f(r,\theta) = (1- r^4)\cos(\theta)\quad\quad g(r,\theta) = 0$$

Problema 2. uponga que $a = 2$, $v = 1$ y que las condiciones iniciales son: $$ f(r,\theta) = 0 \quad\quad g(r,\theta) = r^2(2- r)\sin^8\left(\frac{\theta}{2}\right)$$

Fin Modulo 1



In [ ]:



In [ ]:



In [ ]:

Versión antigua (Lo único que tal vez sirva de algo, son los programas)

Tal vez sirvan de algo para sus tareas.

Entonces, primero veamos algunos modos normales del sistema. Por ejemplo(demasiado simplificado),

$$u(r,\theta, t)_{nk} = J_{n}(\lambda_{nk} r)\,\cos(n\theta)\,\cos(\lambda_{nk} v t)$$

La siguiente función se aplica caso simplificado.



In [3]:

    
def tambor(n, r_max, v, kth_zero, nt, t):
    r = np.r_[0:r_max:100j]
    angle = np.r_[0:2*np.pi:200j]
    ceros = special.jn_zeros(0, nt)
    lamb = ceros[kth_zero]
    u = np.array([special.jn(n, lamb* var_r) * np.cos(n * angle) 
                 * np.cos(lamb * v * t) for var_r in r])
    x = np.array([var_r * np.cos(angle) for var_r in r])
    y = np.array([var_r * np.sin(angle) for var_r in r])
    return x, y, u

Entonces, por ejemplo si $n = 1$, $a = 1$, $v = 1$, $k = 1$ y $t= 0$. Este sería el modo de vibración $(n,k)\rightarrow (1,1)$.



In [4]:

    
x, y, u = tambor(1, 1, 1, 0, 15, 0)



In [5]:

    
plt.figure(figsize = (6, 5))
plt.pcolor(x, y, u, cmap = 'viridis')
plt.axis('off')
plt.colorbar()
plt.show()

Ahora, veamos como lucen todos demás modos de vibración $(n,k)$.



In [6]:

    
def tambor_nk(t = 0, n = 0, kth=0):
    fig = plt.figure(figsize = (6,5))
    ax = fig.add_subplot(1, 1, 1)
    x, y, u = tambor(n, 1, 1, kth, 15, t)
    im = ax.pcolor(x, y, u, cmap = 'viridis')
    ax.set_xlim(xmin = -1.1, xmax = 1.1)
    ax.set_ylim(ymin = -1.1, ymax = 1.1)
    plt.colorbar(im)
    plt.axis('off')
    fig.canvas.draw()
interact_manual(tambor_nk, t = (0, 15,.01), n = (0, 10, 1), kth = (0, 10, 1));

Ahora, tal vez nos interesaría conocer el comportamiento de la membrana cuando sumamos sobre un conjunto de modos $k$. Es decir,

$$u(r,\theta, t)_{n} =\sum_{k = 1}u(r,\theta, t)_{nk} = \sum_{k = 1}J_{n}(\lambda_{nk} r)\,\cos(n\theta)\,\cos(\lambda_{nk} v t) $$

La manera usual de hacer esto es considerar la suma en series de Fourier, es decir a esta suma le falta un coeficiente $A_{nk}$, pero por simplicidad aquí no vamos a considerar este término.

Una posible función para realizar esto sería,



In [7]:

    
def tambor_n_allk(n, r_max, v, nk_zeros, t):
    r = np.r_[0:r_max:100j]
    angle = np.r_[0:2*np.pi:200j]
    ceros = special.jn_zeros(0, nk_zeros)
    lamb = ceros[0]
    u0 = np.array([special.jn(n, lamb* var_r) * np.cos(n * angle)
                  * np.cos(lamb * v * t) for var_r in r])
    for cero in range(1, nk_zeros):
        lamb = ceros[cero]
        u = np.array([special.jn(n, lamb* var_r) * np.cos(n * angle) 
                 * np.cos(lamb * v * t) for var_r in r])
        u0 += u
    x = np.array([var_r * np.cos(angle) for var_r in r])
    y = np.array([var_r * np.sin(angle) for var_r in r])
    return x, y, u0



In [8]:

    
def tambor_n(t = 0, n = 0):
    fig = plt.figure(figsize = (6,5))
    ax = fig.add_subplot(1, 1, 1)
    x, y, u = tambor_n_allk(n, 1, 1, 15, t)
    im = ax.pcolor(x, y, u, cmap = 'viridis')
    ax.set_xlim(xmin = -1.1, xmax = 1.1)
    ax.set_ylim(ymin = -1.1, ymax = 1.1)
    plt.colorbar(im)
    plt.axis('off')
    fig.canvas.draw()
interact_manual(tambor_n, t = (0, 15,.01), n = (0, 10, 1));

Por último, nos queda el caso cuando sumamos sobre todos los modos $n$. Es decir,



In [ ]:



In [37]:

    
def order_n(n, ceros, nk_zeros, angle, v, r, t):
    lamb = ceros[0]
    u0 = np.array([special.jn(n, lamb* var_r) * np.cos(n * angle)
                      * np.cos(lamb * v * t) for var_r in r])
    for cero in range(1, nk_zeros):
        lamb = ceros[cero]
        u = np.array([special.jn(n, lamb* var_r) * np.cos(n * angle) 
                     * np.cos(lamb * v * t) for var_r in r])
        u0 += u
    return u0



In [38]:

    
def tambor(orden_n, r_max, v, nk_zeros, t):
    r = np.r_[0:r_max:100j]
    angle = np.r_[0:2*np.pi:100j]
    ceros = special.jn_zeros(0, nk_zeros)
    u0 = order_n(0, ceros, nk_zeros,  angle, v, r, t)
    for n in range(1, orden_n):
        u = order_n(n, ceros, nk_zeros, angle, v, r, t)
        u0 += u
    x = np.array([var_r * np.cos(angle) for var_r in r])
    y = np.array([var_r * np.sin(angle) for var_r in r])
    return x, y, u0



In [63]:

    
x, y, u = tambor(10, 1, 1, 5, 5)



In [64]:

    
plt.figure(figsize = (5, 5))
plt.pcolor(x, y, u, cmap = 'inferno')
plt.axis('on')
plt.show()