title: fast method with numpy for 2D Heat equation
body: <p>I'm looking for a method for solve the 2D heat equation with python. I have already implemented the finite difference method but is slow motion (to make 100,000 simulations takes 30 minutes). The idea is to create a code in which the end can write,</p>
  
  <pre><code>for t in TIME:  
      DeltaU=f(U)
      U=U+DeltaU*DeltaT
      save(U)
  </code></pre>
  
  <p>How can I do that?</p>
  
  <p>In the first form of my code, I used the 2D method of finite difference, my grill is 5000x250 (x, y). Now I would like to decrease the speed of computing and the idea is to find</p>
  
  <p>DeltaU = f(u)</p>
  
  <p>where U is a heat function. For implementation I used this source <a href="http://www.timteatro.net/2010/10/29/performance-python-solving-the-2d-diffusion-equation-with-numpy/" rel="nofollow">http://www.timteatro.net/2010/10/29/performance-python-solving-the-2d-diffusion-equation-with-numpy/</a> for 2D case, but the run time is more expensive for my necessity. Are there some methods to do this?</p>
  
  <p>Maybe I must to work with the matrix </p>
  
  <pre><code>A=1/dx^2 (2 -1  0  0 ... 0 
           -1  2 -1  0 ... 0
           0  -1  2 -1 ... 0
           .            .
           .            .
           .            .
           0 ...        -1 2)
  </code></pre>
  
  <p>but how to make this in the 2D problem? How to inserting Boundary conditions in A?
  This is the code for the finite difference that I used:</p>
  
  <pre><code>Lx=5000   # physical length x vector in micron
  Ly=250   # physical length y vector in micron
  Nx = 100     # number of point of mesh along x direction
  Ny = 50     # number of point of mesh along y direction
  a = 0.001 # diffusion coefficent
  dx = 1/Nx
  dy = 1/Ny
  dt = (dx**2*dy**2)/(2*a*(dx**2 + dy**2)) # it is 0.04
  x = linspace(0.1,Lx, Nx)[np.newaxis] # vector to create mesh
  y = linspace(0.1,Ly, Ny)[np.newaxis] # vector to create mesh
  I=sqrt(x*y.T) #initial data for heat equation
  u=np.ones(([Nx,Ny])) # u is the matrix referred to heat function
  steps=100000
  for m in range (0,steps):
          du=np.zeros(([Nx,Ny]))
  
          for i in range (1,Nx-1):
  
              for j in range(1,Ny-1):
  
                 dux = ( u[i+1,j] - 2*u[i,j] + u[i-1, j] ) / dx**2
                 duy = ( u[i,j+1] - 2*u[i,j] + u[i, j-1] ) / dy**2            
                 du[i,j] = dt*a*(dux+duy)
  
  
  
  
      # Boundary Conditions
      t1=(u[:,0]+u[:,1])/2
      u[:,0]=t1
      u[:,1]=t1
      t2=(u[0,:]+u[1,:])/2
      u[0,:]=t2
      u[1,:]=t2
      t3=(u[-1,:]+u[-2,:])/2
      u[-1,:]=t3
      u[-2,:]=t3
      u[:,-1]=1
  
  
      filename1='data_{:08d}.txt'
  
      if m%100==0:
          np.savetxt(filename1.format(m),u,delimiter='\t' )
  </code></pre>
  
  <p>For elaborate 100000 steps the run time is about 30 minutes. I would to optimize this code (with the idea presented in the initial lines) to have a run time about 5/10 minutes or minus. How can I do it?</p>