In [3]:
import numpy as np
import matplotlib as mpl; import matplotlib.pyplot as plt
%matplotlib inline
import os; import pickle
$$\frac{\partial^2 \psi}{\partial{x}^2} + \frac{\partial^2 \psi}{\partial{y}^2} - 2\eta \frac{\partial\psi}{\partial{t}} = \frac{1}{v^2}\frac{\partial^2 \psi}{\partial{t}^2}$$

Utilizando las aproximaciones por diferencias finitas para las segundas derivadas y derivadas centrales para la primera derivada con respecto al tiempo se tiene que:

$$\frac{1}{\Delta{x}^2}\left( \psi_{i+1,j}^n + \psi_{i-1,j}^n + \psi_{i,j+1}^n + \psi_{i,j-1}^n - 4\psi_{i,j}^n \right) - \frac{\eta}{\Delta t} \left(\psi_{i,j}^{n+1} - \psi_{i,j}^{n-1} \right) = \frac{1}{v^2 \Delta t^2}\left( \psi_{i,j}^{n+1} + \psi_{i,j}^{n-1} - 2\psi_{i,j}^{n} \right)$$

De donde se puede resolver para $\psi_{i,j}^{n+1}$ en términos de los valores en $n$ y $n-1$:

$$ \psi_{i,j}^{n+1} = \frac{(2-4r^2)}{1+\beta}\psi_{i,j}^n + \frac{r^2}{1+\beta}\left( \psi_{i+1,j}^n + \psi_{i-1,j}^n + \psi_{i,j+1}^n + \psi_{i,j-1}^n \right) + \frac{\beta-1}{1 + \beta}\psi_{i,j}^{n-1}$$

Siendo:

$$r=v\frac{\Delta t}{\Delta x} \textrm{ y } \beta=\eta\Delta{t}v^2$$

La primera iteración debe tratarse como caso especial y usarse para ella la condición sobre la velocidad inicial. Si se asume que la velocidad inicial es igual a cero entonces se tiene que $\psi_{i,j}^{1}=\psi_{i,j}^{-1}$ y con esto se puede encontrar que:

$$\psi_{i,j}^{1} = \frac{1+\beta}{2} \left( \frac{2-4r^2}{1+\beta} \psi_{i,j}^0 + \frac{r^2}{1+\beta} \left( \psi_{i+1,j}^0 + \psi_{i-1,j}^0 + \psi_{i,j+1}^0 + \psi_{i,j-1}^0 \right) \right)$$

$\eta$ define junto con $v$ la escala de tiempo en la que una onda se desvanece siendo ella igual a $\frac{1}{\eta v^2}$.


In [4]:
L=1.
xmin=-L/2  # en metros
xmax=L/2   # en metros
ymin=-L/2
ymax=L/2
v=2. # m/s
tmax=10.*L/v # aproximadamente 5 ires y venires
Nl=80
Nt=1500.
dt=tmax/Nt
print "dt=",dt
tau=L/v  # tiempo de relajación tal que al recorrer medio L ya se haya atenuado significativamente
print "tau=",tau
dx=(xmax-xmin)/Nl
dy=dx
eta=1./v**2/tau
r=v*dt/dx
beta=eta*v**2*dt
print "r=",r
print "beta=",beta
stdev=2*dx # tamaño de las "gotas"


dt= 0.00333333333333
tau= 0.5
r= 0.533333333333
beta= 0.00666666666667

In [5]:
def drop(amp):
    global xmesh
    global ymesn
    global stdev
    x0=np.random.uniform(xmin,xmax)
    y0=np.random.uniform(ymin,ymax)
    #x0=0.
    #y0=0.
    return amp*np.exp(-(xmesh-x0)**2/(2*stdev**2))*np.exp(-(ymesh-y0)**2/(2*stdev**2))
$$ \psi_{i,j}^{n+1} = \frac{(2-4r^2)}{1+\beta}\psi_{i,j}^n + \frac{r^2}{1+\beta}\left( \psi_{i+1,j}^n + \psi_{i-1,j}^n + \psi_{i,j+1}^n + \psi_{i,j-1}^n \right) + \frac{\beta-1}{1 + \beta}\psi_{i,j}^{n-1}$$

In [6]:
def onestep():
    global drum
    drum[0], drum[1] = drum[1], (2.-4.*r**2)/(1.+beta)*drum[1]-\
    (1.-beta)/(1.+beta)*drum[0]+\
    r**2/(1.+beta)*(np.roll(drum[1],1,axis=1)+\
                    np.roll(drum[1],-1,axis=1)+\
                    np.roll(drum[1],1,axis=0)+\
                    np.roll(drum[1],-1,axis=0))
def firststep():
    global drum
    global eta
    global dt
    global r
    drum[1] = (2.-4.*r**2)/(1.+beta)*drum[0]+r**2/(1.+beta)*(np.roll(drum[0],1,axis=0)+np.roll(drum[0],-1,axis=0)+np.roll(drum[0],1,axis=1)+np.roll(drum[0],-1,axis=1))
def fixborders():
    global drum
    drum[1,0]=0
    drum[1,-1]=0
    drum[1,:,-1]=0
    drum[1,:,0]=0

In [7]:
#extremos fijos
xcoords=np.linspace(xmin,xmax,Nl)
ycoords=np.linspace(ymin,ymax,Nl)
xmesh,ymesh=np.meshgrid(xcoords,ycoords)
drum=np.zeros((2,Nl,Nl))
drum[0]=drop(0.01)
fixborders()
firststep()
fixborders()

In [8]:
query=[]

In [9]:
# Animación
%matplotlib osx
from matplotlib import animation

fig=plt.figure(figsize=(10,10))
plt.hold(True)
tplot=plt.imshow(drum[0],cmap='winter',vmin=0.,vmax=0.02,interpolation='None')
#plt.xlim(xmin+dx,xmax-dx)
#plt.ylim(ymin+dy,ymax-dy)

def animate(i):
    global drum
    global query
    onestep()
    query.append(drum[1,Nl/2,Nl/2])
    if np.random.random()<=10*100./Nt:
        drum[1]+=drop(0.01*np.random.random())
    fixborders()
    onestep();fixborders()
    onestep();fixborders()
    tplot.set_array(np.abs(drum[0]))
    return 0

animacion = animation.FuncAnimation(fig, animate, np.arange(2,1e6),interval=1, blit=False)
plt.show()
plt.hold(False)

In [75]:
len(query)


Out[75]:
4985

In [79]:
%matplotlib inline
plt.figure()
plt.hist(query,bins=30,normed=True)
plt.show()



In [25]:
%%timeit
onestep()


10000 loops, best of 3: 142 µs per loop

In [19]:
#Movie
plt.ioff()
pad="000"
for i in range(1,int(Nt)):
    print i,
    plt.figure(figsize=(10,10))
    onestep()
    if np.random.randint(1,10)>8:
        drum[1]+=drop(0.5*np.random.random())
    fixborders()
    onestep();fixborders()
    onestep();fixborders()
    tplot=plt.pcolor(xcoords,ycoords,np.abs(drum[0]),cmap='gist_heat',vmin=0.,vmax=1.)
    plt.xlim(xmin+dx,xmax-dx)
    plt.ylim(ymin+dy,ymax-dy)
    if i>=10:
        pad="00"
    if i>=100:
        pad="0"
    if i>=1000:
        pad=""
    tplot.axes.get_xaxis().set_visible(False)
    tplot.axes.get_yaxis().set_visible(False)
    plt.savefig("/Users/juan/Desktop/pool/water-"+pad+str(i)+".png")
    plt.close()


1 2 3 4 5 6 7 8
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-19-e42e53a355ba> in <module>()
     22     tplot.axes.get_xaxis().set_visible(False)
     23     tplot.axes.get_yaxis().set_visible(False)
---> 24     plt.savefig("/Users/juan/Desktop/pool/water-"+pad+str(i)+".png")
     25     plt.close()

/Users/juan/anaconda/lib/python2.7/site-packages/matplotlib/pyplot.pyc in savefig(*args, **kwargs)
    575 def savefig(*args, **kwargs):
    576     fig = gcf()
--> 577     res = fig.savefig(*args, **kwargs)
    578     draw()   # need this if 'transparent=True' to reset colors
    579     return res

/Users/juan/anaconda/lib/python2.7/site-packages/matplotlib/figure.pyc in savefig(self, *args, **kwargs)
   1474             self.set_frameon(frameon)
   1475 
-> 1476         self.canvas.print_figure(*args, **kwargs)
   1477 
   1478         if frameon:

/Users/juan/anaconda/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2209                 orientation=orientation,
   2210                 bbox_inches_restore=_bbox_inches_restore,
-> 2211                 **kwargs)
   2212         finally:
   2213             if bbox_inches and restore_bbox:

/Users/juan/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    531             _png.write_png(renderer._renderer.buffer_rgba(),
    532                            renderer.width, renderer.height,
--> 533                            filename_or_obj, self.figure.dpi)
    534         finally:
    535             if close:

KeyboardInterrupt: 

In [ ]:


In [ ]: