In [3]:
import numpy as np
import matplotlib as mpl; import matplotlib.pyplot as plt
%matplotlib inline
import os; import pickle
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"
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))
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]:
In [79]:
%matplotlib inline
plt.figure()
plt.hist(query,bins=30,normed=True)
plt.show()
In [25]:
%%timeit
onestep()
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()
In [ ]:
In [ ]: