notebook.community

Edit and run

Heat Equation

\begin{eqnarray*} \frac{\partial^2 u(x,y,t)}{\partial t^2} & = & c^2\nabla^2 u(x,y,t) =c^2\left( \ \frac{\partial^2 u(x,y,t)}{\partial x^2}+\frac{\partial^2 u(x,y,t)}{\partial y^2}\right) \end{eqnarray*}

with initial conditions: $$ u(x,y,0)=f(x,y), u_t(x,y,0)=g(x,y)$$ and homogeneous Neumann boundary condition: $$\vec n\cdot \nabla u|_{\partial \Omega}=0$$

Numeric Scheme

  • discretization: $u(t_m,x_i,y_j)= U^m_{i,j}$

  • Equations \begin{eqnarray} \frac{\partial^2 u}{\partial t^2} &\sim& \frac{U^{k+1}_{i,j}-2U^k_{i,j}+U^{k-1}_{i,j}}{\triangle t^2}\\ \frac{\partial^2 u}{\partial x^2} &\sim& \frac{ U^k_{i+1,j}-2U^k_{i,j}+U^k_{i-1,j}}{\triangle x^2}\\ \frac{\partial^2 u}{\partial y^2} &\sim& \frac{ U^k_{i,j+1}-2U^k_{i,j}+U^k{i,j-1}}{\triangle y^2}\\ \frac{U^{k+1}_{i,j}-2U^k_{i,j}+U^{k-1}_{i,j}}{\triangle t^2} &=& c^2\left( \frac{ U^k_{i+1,j}-2U^k_{i,j}+U^k_{i-1,j}}{\triangle x^2}+ \frac{ U^k_{i,j+1}-2U^k_{i,j}+U^k_{i,j-1}}{\triangle y^2} \right) \end{eqnarray}

  • scheme (main steps): \begin{eqnarray} \gamma_x,\gamma_y,\gamma &=& \frac{c\Delta t}{\Delta x}, \frac{c\Delta t}{\Delta y}, 2(1-\gamma_x^2-\gamma_y^2)\\ U_{i,j}^{k+1} &=& \gamma U_{i,j}^k+\gamma_x^2[U^k_{i,j+1}+U^k_{i,j-1}]+\gamma_y^2[U_{i+1,j}^k+U_{i-1,j}^k]-U_{i,j}^{k-1} \end{eqnarray}

  • Initial Condition, $ u(x,y,0)=f(x,y), u_t(x,y,0)=g(x,y)$: \begin{eqnarray} U_{i,j}^0 &=& f_{i,j},\quad \frac{U_{ij}^1-U_{i,j}^{-1}}{2\Delta t}=g_{i,j}\\ \to U_{i,j}^{-1} &=& U_{i,j}^1 - 2\Delta t g_{i,j} \\ \frac{2 U_{i,j}^{1}-2f_{i,j}-2\Delta t g_{i,j}}{\Delta t^2} &=& c^2\left( \frac{f_{i,j+1}-2f_{i,j}+f_{i,j-1}}{\Delta x^2} +\frac{f_{i+1,j}-2f_{i,j}+f_{i-1,j}}{\Delta y^2}\right) \\ U_{i,j}^1 &=& \frac{\gamma}{2}f_{i,j}+\frac{\gamma_x^2}{2}(f_{i,j+1}+f_{i,j-1}) +\frac{\gamma_y^2}{2}(f_{i+1,j} +f_{i-1,j})+\Delta t g_{i,j} \\ \end{eqnarray}

  • Boundary conditions (rectangular region), $\vec n\cdot \nabla u|_{\partial \Omega}=0$:
    • the left boundary, $j=0$ $$U_{i,0}^{k+1} = \gamma U_{i,0}^k+2\gamma_x^2 U^k_{i,j-1}+\gamma_y^2[U_{i+1,0}^k+U_{i-1,0}]-U_{i,0}^{k-1}$$
    • the right boundary, $j=n$ $$ U_{i,n}^{k+1} = \gamma U_{i,n}^k+2\gamma_x^2U^k_{i,n-1}+\gamma_y^2[U_{i+1,n}^k+U_{i-1,n}]-U_{i,n}^{k-1}$$
    • the bottom boundary, $i=0$ $$ U_{0,j}^{k+1} = \gamma U_{0,j}^k+\gamma_x^2[U^k_{0,j+1}+U^k_{0,j-1}]+2\gamma_y^2U_{1,j}^k-U_{0,j}^{k-1}$$
    • the top boundary, $i=m$ $$ U_{m,j}^{k+1} = \gamma U_{m,j}^k+\gamma_x^2[U^k_{m,j+1}+U^k_{m,j-1}]+2\gamma_y^2U_{1,j}^k-U_{m,j}^{k-1}$$
  • Corner points:
    • the upper right corner, $i=m$ and $j=n$ $$U_{m,n}^{k+1} = \gamma U_{m,n}^k+2\gamma_x^2U^k_{m,n-1}+2\gamma_y^2 U_{m-1,n}-U_{m,n}^{k-1} $$
    • the lower right corner, $i=0$ and $j=n$ $$U_{0,n}^{k+1} = \gamma U_{0,n}^k+2\gamma_x^2U^k_{0,n-1}+2\gamma_y^2 U_{1,n}-U_{0,n}^{k-1} $$
    • the upper left corner, $i=m$ and $j=0$ $$U_{m,0}^{k+1} = \gamma U_{m,0}^k+2\gamma_x^2 U^k_{m,1}+2\gamma_y^2U_{m-1,0}^k-U_{m,0}^{k-1} $$
    • the lower left corner, $i=0$ and $j=0$ $$U_{0,0}^{k+1} = \gamma U_{0,0}^k+2\gamma_x^2 U^k_{0,1}+2\gamma_y^2U_{1,0}^k-U_{0,0}^{k-1}$$

Stability

$$\triangle t \le \frac{1}{2c^2}\frac{(\triangle x \triangle y)^2}{\triangle x^2+\triangle y^2}$$ 




In [ ]:



    






    











In [24]:



    


from sympy import Symbol,symbols,solve,simplify,init_printing
from sympy.physics.mechanics import dynamicsymbols



    











In [2]:



    


def WaveScheme():
    ut0,ut1,ut2=symbols("ut0 ut1 ut2")
    ux=symbols('ux:' + str(3))
    uy=symbols('uy:' + str(3))
    dt,dx,dy,u,s=symbols("dt dx dy u s")
    i,j,k=symbols("i j k", type=int)
    EqS=(ut0-2*u+ut2)/dt**2.-((ux[0]-2*u+ux[2])/dx**2.+(uy[0]-2*u+uy[2])/dy**2.)
    result=solve(EqS,ut2)
    result1=result[0].subs({ut0:str('U(i,j,k-1)'),u:str('U(i,j,k)'),ux[0]:str('U(i-1,j,k)'),ux[2]:str('U(i+1,j,k)'), \
                  uy[0]:str('U(i,j-1,k)'),uy[2]:str('U(i,j+1,k)')})
    #result=solve(EqS1,ut2)
    return result1



    











In [22]:



    


def WaveScheme2(c):
    ut0,ut1,ut2=symbols("ut0 ut1 ut2")
    ux=symbols('ux:' + str(3))
    uy=symbols('uy:' + str(3))
    dt,dx,dy,u,s=symbols("dt dx dy u s")
    i,j,k=symbols("i j k")
    
    dudt2=(ut0-2*u+ut2)/dt**2.
    dudx2=(ux[0]-2*u+ux[2])/dx**2.+(uy[0]-2*u+uy[2])/dy**2.
    

    EqS=dudt2-c**2*dudx2
    result=solve(EqS,ut2)
    result1=result[0].subs({ut0:str('U(i,j,k-1)'),u:str('U(i,j,k)'),ux[0]:str('U(i-1,j,k)'),ux[2]:str('U(i+1,j,k)'), \
                  uy[0]:str('U(i,j-1,k)'),uy[2]:str('U(i,j+1,k)')})
    #result=solve(EqS1,ut2)
    return result1



    











In [18]:



    


init_printing(True)
scheme1=WaveScheme();scheme1



    


















    
Out[18]:







$$\frac{1}{dx^{2} dy^{2}} \left(dt^{2} dx^{2} \left(- 2.0 U{\left (i,j,k \right )} + U{\left (i,j - 1,k \right )} + U{\left (i,j + 1,k \right )}\right) + dt^{2} dy^{2} \left(- 2.0 U{\left (i,j,k \right )} + U{\left (i - 1,j,k \right )} + U{\left (i + 1,j,k \right )}\right) + dx^{2} dy^{2} \left(2.0 U{\left (i,j,k \right )} - U{\left (i,j,k - 1 \right )}\right)\right)$$
















In [25]:



    


init_printing(True)
c=Symbol("c")
scheme1=WaveScheme2(c);scheme1



    


















    
Out[25]:







$$\frac{1}{dx^{2} dy^{2}} \left(c^{2} dt^{2} dx^{2} \left(- 2.0 U{\left (i,j,k \right )} + U{\left (i,j - 1,k \right )} + U{\left (i,j + 1,k \right )}\right) + c^{2} dt^{2} dy^{2} \left(- 2.0 U{\left (i,j,k \right )} + U{\left (i - 1,j,k \right )} + U{\left (i + 1,j,k \right )}\right) + dx^{2} dy^{2} \left(2.0 U{\left (i,j,k \right )} - U{\left (i,j,k - 1 \right )}\right)\right)$$
















In [21]:



    






    











In [ ]:



    






    











In [3]:



    


%pylab inline
import numpy as np
import matplotlib.pyplot as plt
from numpy import exp,sin,pi
from matplotlib import animation
from JSAnimation import IPython_display
import time



    


















    






Populating the interactive namespace from numpy and matplotlib

















In [2]:



    


def MainScheme(uk,ukm1,dt,dx,c):
    gx2 = c*dt/dx
    gy2 = c*dt/dy
    gamma = 2*(1 - gx2 - gy2)
    u[:,:,1] = 0.5*gamma*uk-dt-ulm1
    u[1:-1,1:-1,1] += 0.5*gx2*(uk[1:-1,2:]+uk[1:-1,:-2])
    u[1:-1,1:-1,1] += 0.5*gy2*(uk[:-2,1:-1]+uk[2:,1:-1])
    return u



    











In [2]:



    


class wave2d(object):
    def __init__(self,height,width,T,nx,ny,nt,c):
        self.x = np.linspace(-0.5*width,0.5*width,nx+1)
        self.y = np.linspace(-0.5*height,0.5*height,ny+1)
        self.t = np.linspace(0,T,nt+1)

        self.dx = self.x[1]-self.x[0]
        self.dy = self.y[1]-self.y[0]
        self.dt = self.t[1]-self.t[0]
        
        self.xx,self.yy = np.meshgrid(self.x,self.y)

        # Gamma_x squared
        self.gx2 = c*self.dt/self.dx

        # Gamma_y squared
        self.gy2 = c*self.dt/self.dy

        # 2*(1-gamma_x^2-gamma_y^2)
        self.gamma = 2*(1 - self.gx2 - self.gy2)
        

    def solve(self,ffun,gfun):

        #f = ffun(self.xx,self.yy)        
        f = ffun(self.xx,self.yy)
        g = gfun(self.xx,self.yy) 

        u = np.zeros((ny+1,nx+1,nt+1))

        # Set initial condition
        u[:,:,0] = f
        U= np.copy( u[:,:,0])

        """ Compute first time step """
        u[:,:,1] = 0.5*self.gamma*f+self.dt*g
        u[1:-1,1:-1,1] += 0.5*self.gx2*(f[1:-1,2:]+f[1:-1,:-2])
        u[1:-1,1:-1,1] += 0.5*self.gy2*(f[:-2,1:-1]+f[2:,1:-1])

        for k in range(1,nt):
            # Every point contains these terms
            u[:,:,k+1] = self.gamma*u[:,:,k] - u[:,:,k-1]
            
            # Interior points
            u[1:-1,1:-1,k+1] += self.gx2*(u[1:-1,2:,k]+u[1:-1,:-2,k]) + \
                                self.gy2*(u[2:,1:-1,k]+u[:-2,1:-1,k]) 

            # Top boundary
            u[0,1:-1,k+1] +=  2*self.gy2*u[1,1:-1,k] + \
                                self.gx2*(u[0,2:,k]+u[0,:-2,k])
                                      
            # Right boundary
            u[1:-1,-1,k+1] += 2*self.gx2*u[1:-1,-2,k] + \
                                self.gy2*(u[2:,-1,k]+u[:-2,-1,k])

            # Bottom boundary
            u[-1,1:-1,k+1] +=  2*self.gy2*u[-2,1:-1,k] + \
                                 self.gx2*(u[-1,2:,k]+u[-1,:-2,k]) 
                               
            # Left boundary
            u[1:-1,0,k+1] += 2*self.gx2*u[1:-1,1,k] + \
                               self.gy2*(u[2:,0,k]+u[:-2,0,k])

            # Top right corner
            u[0,-1,k+1] += 2*self.gx2*u[0,-2,k] + \
                           2*self.gy2*u[1,-1,k]

            # Bottom right corner
            u[-1,-1,k+1] += 2*self.gx2*u[-1,-2,k] + \
                            2*self.gy2*u[-2,-1,k]

            # Bottom left corner
            u[-1,0,k+1] += 2*self.gx2*u[-1,1,k] + \
                           2*self.gy2*u[-2,0,k]
    
            # Top left corner
            u[0,0,k+1] += 2*self.gx2*u[0,1,k] + \
                          2*self.gy2*u[1,0,k]    
        
        return u



    











In [3]:



    


class wave2dR(object):
    def __init__(self,height,width,T,nx,ny,nt,c):
        self.x = np.linspace(-0.5*width,0.5*width,nx+1)
        self.y = np.linspace(-0.5*height,0.5*height,ny+1)
        self.t = np.linspace(0,T,nt+1)

        self.dx = self.x[1]-self.x[0]
        self.dy = self.y[1]-self.y[0]
        self.dt = self.t[1]-self.t[0]
        
        self.xx,self.yy = np.meshgrid(self.x,self.y)

        # Gamma_x squared
        self.gx2 = c*self.dt/self.dx

        # Gamma_y squared
        self.gy2 = c*self.dt/self.dy

        # 2*(1-gamma_x^2-gamma_y^2)
        self.gamma = 2*(1 - self.gx2 - self.gy2)
        

    def solve(self,ffun,gfun):

        #f = ffun(self.xx,self.yy)        
        f = ffun(self.xx,self.yy)
        g = gfun(self.xx,self.yy) 

        u = np.zeros((ny+1,nx+1,nt+1))

        # Set initial condition
        
        u[:,:,0] = f
        
        """ Compute first time step """
        u[:,:,1] = 0.5*self.gamma*f+self.dt*g
        u[1:-1,1:-1,1] += 0.5*self.gx2*(f[1:-1,2:]+f[1:-1,:-2])
        u[1:-1,1:-1,1] += 0.5*self.gy2*(f[:-2,1:-1]+f[2:,1:-1])
        
        yo,xo = np.ogrid[-ny/2:ny-ny/2+1, -nx/2:nx-nx/2+1]
        mask = xo**2 + yo**2 >= nx*nx/4

        array = u[:,:,1]
        array[mask] = 0.
        u[:,:,1] = array
        print(len(u[:,:,1]))
        for k in range(1,nt):
            # Every point contains these terms
            u[:,:,k+1] = self.gamma*u[:,:,k] - u[:,:,k-1]
            # Interior points
            u[1:-1,1:-1,k+1] += self.gx2*(u[1:-1,2:,k]+u[1:-1,:-2,k]) + \
                                self.gy2*(u[2:,1:-1,k]+u[:-2,1:-1,k]) 
            #UU= u[:,:,k+1]
            #UU[mask]=0.
            #u[:,:,k+1]=UU
            u[mask,k+1]=0.
        return u



    











In [16]:



    


T = 4
height = 2
width = 4
c = 1
nt = 800
nx = 250
ny = 125



    











In [18]:



    


nt=800
wave_eq = wave2d(height,width,T,nx,ny,nt,c)

# Initial value functions
f = lambda x,y: np.exp(-10*(x**2+y**2))
g = lambda x,y: 0 

u = wave_eq.solve(f,g)



    











In [19]:



    


x = wave_eq.x
y = wave_eq.y

plt.imshow(u[:,:,-1],extent=[x[0],x[-1],y[0],y[-1]])



    


















    
Out[19]:








<matplotlib.image.AxesImage at 0x11d7e13d0>










    
























In [28]:



    


dx=width/float(nx)
dy=height/float(ny)

sqrt(1/c**2*dx**2*dy**2/(dx**2+dy**2)/2.)



    


















    
Out[28]:








0.0080000000000000002

















In [26]:



    


4/800.



    


















    
Out[26]:








0.005

Question

What about the graphical output if we set the nt=400 or less?





In [ ]:



    






    











In [20]:



    


anim = plt.figure(figsize=(16,8))
ax = anim.add_subplot(111)
#ax = anim.add_subplot(111)
#ax.set_title("2D Wave Equation",fontsize=14)
T = 4
height = 2
width = 4
c = 1
nt = 800
nx = 250
ny = 125

wave_eq = wave2d(height,width,T,nx,ny,nt,c)

# Initial value functions
f = lambda x,y: np.exp(-10*(x**2+y**2))
g = lambda x,y: 0 

u = wave_eq.solve(f,g)

x = wave_eq.x
y = wave_eq.y

def init():
    return ax.imshow(u[:,:,0],vmin=0,vmax=1,extent=[x[0],x[-1],y[0],y[-1]])

def animate(i):
    return ax.imshow(u[:,:,10*i],vmin=0,vmax=1,extent=[x[0],x[-1],y[0],y[-1]])

 
animation.FuncAnimation(anim, animate, init_func=init, frames=nt/10, interval=5)



    


















    
Out[20]:











    
    

    
    

    
    
    
    
    
    
    
    
    
  

     Once 
     Loop 
     Reflect 
  























In [167]:

Drum Wave





In [ ]:



    






    











In [10]:



    


T = 1
height = 4
width = 4
c = 1
nt = 200
nx = 250
ny = 250



    











In [11]:



    


#wave_eq = wave2dR(height,width,T,nx,ny,nt,c)

# Initial value functions
#f = lambda x,y:  np.exp(-x**2.-y**2.)-1. if x**2.+y**2.<=4. else 0;
#"""
"""
def f(x,y):    
    U=0*x
    if (np.any((x)**2+(y)**2)<=4):
       U=np.exp(-(x)**2.-(y)**2.)-exp(-4)
    return U 
"""
def f(x,y):    
    R=x**2+y**2
    U=np.cos(np.pi*np.sqrt((x)**2.+(y)**2.)/4.)
    U[R>=2.]=0
    return U 
g = lambda x,y: 0 
#def g(x,y): return 0.
wave_eq = wave2dR(height,width,T,nx,ny,nt,c)

u = wave_eq.solve(f,g)



    


















    






251

















In [17]:



    


x = wave_eq.x
y = wave_eq.y

plt.imshow(u[:,:,50],extent=[x[0],x[-1],y[0],y[-1]])
plt.colorbar()



    


















    
Out[17]:








<matplotlib.colorbar.Colorbar instance at 0x123444560>










    
























In [49]:



    


13952/251



    


















    
Out[49]:








55

















In [3]:



    


xx=np.linspace(-0.5*width,0.5*width,nx+1)
yy=np.linspace(-0.5*width,0.5*width,nx+1)

xm,ym = np.meshgrid(xx,yy)
mask= xm**2+ym**2
mask=np.meshgrid(xx,yy)



    











In [6]:



    


n=7
r=3
y,x = np.ogrid[-r:n-r, -r:n-r]
mask = x*x + y*y <= r*r
arr=array*5
array = np.ones((n, n))
array[mask] = 255
array



    


















    
Out[6]:








array([[   1.,    1.,    1.,  255.,    1.,    1.,    1.],
       [   1.,  255.,  255.,  255.,  255.,  255.,    1.],
       [   1.,  255.,  255.,  255.,  255.,  255.,    1.],
       [ 255.,  255.,  255.,  255.,  255.,  255.,  255.],
       [   1.,  255.,  255.,  255.,  255.,  255.,    1.],
       [   1.,  255.,  255.,  255.,  255.,  255.,    1.],
       [   1.,    1.,    1.,  255.,    1.,    1.,    1.]])

















In [23]:



    


nx=12
r=3
y,x = np.ogrid[-nx/2:nx-nx/2+1, -nx/2:nx-nx/2+1]
mask = x*x + y*y >= nx/2*nx/2

array = np.ones((nx+1, nx+1))
array[mask] = 0
array



    


















    
Out[23]:








array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

















In [24]:



    


x*x



    


















    
Out[24]:








array([[-6, -5, -4, -3, -2, -1,  0,  1,  2,  3,  4,  5,  6]])

















In [9]:



    


anim = plt.figure(figsize=(16,8))
ax = anim.add_subplot(111)
#ax = anim.add_subplot(111)
#ax.set_title("2D Wave Equation",fontsize=14)
T=8
nt=int(T/0.005)
height = 4
width = 4
c = 1

nx = 400
ny = 400
wave_eq = wave2dR(height,width,T,nx,ny,nt,c)

# Initial value functions
def f(x,y):    
    R=x**2+y**2
    U=np.cos(np.pi*np.sqrt((x)**2.+(y)**2.)/4.)
    U[R>=2.]=0
    return U  
    
g = lambda x,y: 0 
#def g(x,y): return 0.

u = wave_eq.solve(f,g)

x = wave_eq.x
y = wave_eq.y
frate=nt/40

def init():
    #return ax.imshow(u[:,:,0],vmin=-2,vmax=2,extent=[x[0],x[-1],y[0],y[-1]])
    return ax.imshow(u[:,:,0],vmin=-1,vmax=1,extent=[x[0],x[-1],y[0],y[-1]])
def animate(i):
    return ax.imshow(u[:,:,i*frate],vmin=-1,vmax=1,extent=[x[0],x[-1],y[0],y[-1]])

 
animation.FuncAnimation(anim, animate, init_func=init, frames=nt/frate, interval=100)



    


















    






401










    
Out[9]:











    
    

    
    

    
    
    
    
    
    
    
    
    
  

     Once 
     Loop 
     Reflect 
  























In [ ]:



    






    











In [1]:



    


from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm

import numpy as np
from numpy import sin, cos, arctan, arctan2, array, sqrt, pi, linspace, meshgrid
import scipy
import scipy.special
from scipy.special import jn, jn_zeros
import time, sys, os

makeplots = True



    











In [2]:



    


v = 1
fns = 'cc'
ms = (1,2)
#ms = (1,2,3)
ns = (0,1)
#ns = (0,1,2)
    
f1 = {'s':sin,'c':cos}[fns[0]]
f2 = {'s':sin,'c':cos}[fns[1]]



    











In [3]:



    


def k(m,n):
    return jn_zeros(n,m)[m-1] # m is 0-indexed here


def generate(X, Y, t, m, n, v):
    theta = arctan2(Y,X)
    R = sqrt(X**2 + Y**2)
    # We know z = J_n(k*r)*      cos(n*theta)*cos(k*v*t)
    result =      jn(n,k(m,n)*R)*f1(n*theta)* f2(k(m,n)*v*t)
    result[R>1] = 0  # we plot points from the square, but physically require this.
    return result



    











In [4]:



    


#plt.ion()
fig = plt.figure()
# It seems to me that things look best if you have 6 inches in width
# for each m and 2 in in height for each n.

fig.set_size_inches([6*len(ms),2*len(ns)])
axs = {}
rows, cols = len(ns), 2*len(ms)



    











In [5]:



    


idx = 1
for m in ms:
    axs[m] = {}
    for n in ns:
        axs[m][n] = (fig.add_subplot(rows,cols,idx, projection='3d'),
                    fig.add_subplot(rows,cols,idx+1))
        idx += 2
if 0:
    axs[1] = {0: (fig.add_subplot(rows,cols,1, projection='3d'), fig.add_subplot(rows,cols,2)),
              1: (fig.add_subplot(rows,cols,3, projection='3d'), fig.add_subplot(rows,cols,4)),}
    axs[2] = {0: (fig.add_subplot(rows,cols,5, projection='3d'), fig.add_subplot(rows,cols,6)),
              1: (fig.add_subplot(rows,cols,7, projection='3d'), fig.add_subplot(rows,cols,8)),}



    











In [6]:



    


xs = linspace(-1, 1, 100)
ys = linspace(-1, 1, 100)
X, Y = meshgrid(xs, ys)
#Z = generate(X, Y, 0.0)



    











In [7]:



    


wframes = {}
contours = {}
for m in ms:
    wframes[m] = {}
    contours[m] = {}
    for n in ns:
        wframes[m][n] = None
        contours[m][n] = None

if makeplots:
    dname = fns[0]+fns[1] + '.'.join([str(i) for i in ms]) + '_' + '.'.join([str(i) for i in ns])
    if not os.path.exists(dname):
        os.makedirs(dname)
tstart = time.time()



    











In [8]:



    


first = True
frames_per = 50
periods = 2
# Note: unlike sin and cos, Jn's zeros are not integer multiples of
# each other.  Therefore, this loop goes over a defined number of
# periods of the lowest mode, but others won't fit evenly.
junk = raw_input('go for barney')
for (idx,t) in enumerate(linspace(0, periods*2*pi/jn_zeros(0,1)[0], periods*frames_per)):
    if not first:
        for m in ms:
            for n in ns:
                axs[m][n][0].collections.remove(wframes[m][n])
                for c in contours[m][n].collections:
                    axs[m][n][1].collections.remove(c)
                    
    first = False

    for m in ms:
        for n in ns:
            Z = generate(X, Y, t, m, n, v)
            wframes[m][n] = axs[m][n][0].plot_surface(X, Y, Z, rstride=4, cstride=4, alpha=0.3)

            if n == 0:
                axs[m][n][0].set_zlim(-0.7,.7)
                axs[m][n][1].imshow(Z,vmin=-0.7,vmax=0.7)
            else:
                axs[m][n][0].set_zlim(-0.5,0.5)
                axs[m][n][1].imshow(Z,vmin=-0.5,vmax=0.5)
            # The funny business with levels here is because you won't
            # get a contour exactly at zero that necessarily tracks
            # around both sides of the circle due to the fact that
            # we've discretized things.
            levels = [-0.000000001,0.0,0.000000001]
            contours[m][n] = axs[m][n][1].contour(Z, levels, colors='k',
                                                  linestyles='solid', linewidths=2)
    plt.draw()
    if makeplots:
        plt.savefig( 'drum/img' + '%04d'%(idx) + '.png')

#print ('FPS: %f' % (periods*frames_per / (time.time() - tstart)))



    


















    






go for barneypos

















In [29]:



    


%%bash

ffmpeg     -i drum/img%04d.png  movie.webm
#ffmpeg -q:a 5 -r 5 -b:v 19200 -i drum/img%04d.png movie.mp4
#ffmpeg -h



    


















    






ffmpeg version 2.2.4-tessus Copyright (c) 2000-2014 the FFmpeg developers
  built on Jun 29 2014 16:35:46 with clang version 3.3 (tags/RELEASE_33/final)
  configuration: --cc=/opt/local/bin/clang-mp-3.3 --prefix=/Users/tessus/data/ext/ffmpeg/sw --as=yasm --extra-version=tessus --disable-shared --enable-static --disable-ffplay --enable-gpl --enable-pthreads --enable-postproc --enable-libmp3lame --enable-libtheora --enable-libvorbis --enable-libx264 --enable-libx265 --enable-libxvid --enable-libspeex --enable-bzlib --enable-zlib --enable-libopencore-amrnb --enable-libopencore-amrwb --enable-libxavs --enable-version3 --enable-libvo-aacenc --enable-libvo-amrwbenc --enable-libvpx --enable-libgsm --enable-libopus --enable-libmodplug --enable-fontconfig --enable-libfreetype --enable-libass --enable-libbluray --enable-filters --disable-indev=qtkit --enable-runtime-cpudetect
  libavutil      52. 66.100 / 52. 66.100
  libavcodec     55. 52.102 / 55. 52.102
  libavformat    55. 33.100 / 55. 33.100
  libavdevice    55. 10.100 / 55. 10.100
  libavfilter     4.  2.100 /  4.  2.100
  libswscale      2.  5.102 /  2.  5.102
  libswresample   0. 18.100 /  0. 18.100
  libpostproc    52.  3.100 / 52.  3.100
Input #0, image2, from 'drum/img%04d.png':
  Duration: 00:00:04.00, start: 0.000000, bitrate: N/A
    Stream #0:0: Video: png, rgba, 1200x400 [SAR 3937:3937 DAR 3:1], 25 fps, 25 tbr, 25 tbn, 25 tbc
File 'movie.webm' already exists. Overwrite ? [y/N] Not overwriting - exiting
Conversion failed!

<video src="movie.webm" autoplay,control> Movie </video>





In [ ]:

Content source: cchuang2009/DyLabIPy

Similar notebooks:

notebook.community | gallery | about