title: Finite difference method for 3D diffusion/heat equation
body: <p>I'm trying to use finite differences to solve the diffusion equation in 3D. I think I'm having problems with the main loop. In particular the discrete equation is:</p>
  <p><img src="https://i.imgur.com/F1nyx73.jpg" alt="" /></p>
  <p>With Neumann boundary conditions (in just one face as an example):</p>
  <p><img src="https://i.imgur.com/s1yRDIu.jpg" alt="" /></p>
  <p>Now the code:</p>
  <pre><code>import numpy as np
  from matplotlib import pyplot, cm
  from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots
  %matplotlib inline
  
  kx = 15     #Number of points
  ky = 15
  kz = 15
  largx = 90  #Domain length.
  largy = 90
  largz = 90   
  
  dt4 = 1/2 #Time delta (arbitrary for the time).
  dx4 = largx/(kx-1)    #Position deltas.
  dy4 = largy/(ky-1)
  dz4 = largz/(kz-1) 
  
  Tin = 25    #Initial temperature
  kapp = 0.23
  
  Tamb3d = 150 #Ambient temperature
  
  #Heat per unit of area. One for each face.
  qq1=0 / (largx*largz)
  qq2=0 / (largx*largz)
  qq3=0 / (largz*largy)
  qq4=0 / (largz*largy)
  qq5=0 / (largx*largy)
  qq6=1000 / (largx*largy)
  
  
  x4 = np.linspace(0,largx,kx)   
  y4 = np.linspace(0,largy,ky)
  z4 = np.linspace(0,largz,kz)
  
  
  #Defining the function.
  def diff3d(tt):
      
      w2 = np.ones((kx,ky,kz))*Tin         #Temperature array
      wn2 = np.ones((kx,ky,kz))*Tin
      
      for k in range(tt+2):
          wn2 = w2.copy()
          w2[1:-1,1:-1,1:-1] = (wn2[1:-1,1:-1,1:-1] + 
                          kapp*dt4 / dy4**2 *                                       
                          (wn2[1:-1, 2:,1:-1] - 2 * wn2[1:-1, 1:-1,1:-1] + wn2[1:-1, 0:-2,1:-1]) +  
                          kapp*dt4 / dz4**2 *                                       
                          (wn2[1:-1,1:-1,2:] - 2 * wn2[1:-1, 1:-1,1:-1] + wn2[1:-1, 1:-1,0:-2]) +
                          kapp*dt4 / dx4**2 *
                          (wn2[2:,1:-1,1:-1] - 2 * wn2[1:-1, 1:-1,1:-1] + wn2[0:-2, 1:-1,1:-1]))
  
          #Neumann boundary (dx=dy=dz for the time)
          w2[0,:,:] =   w2[0,:,:] + 2*kapp* (dt4/(dx4**2)) * (w2[1,:,:] - w2[0,:,:] - qq1 * dx4/kapp)
          w2[-1,:,:] =  w2[-1,:,:] + 2* kapp*(dt4/(dx4**2)) * (w2[-2,:,:] - w2[-1,:,:] + qq2 * dx4/kapp)
          w2[:,0,:] =   w2[:,0,:] + 2*kapp* (dt4/(dx4**2)) * (w2[:,1,:] - w2[:,0,:] - qq3 * dx4/kapp)
          w2[:,-1,:] =  w2[:,-1,:] + 2*kapp* (dt4/(dx4**2)) * (w2[:,-2,:] - w2[:,-1,:] + qq4 * dx4/kapp)
          w2[:,:,0] =   w2[:,:,0] + 2 *kapp* (dt4/(dx4**2)) * (w2[:,:,-1] - w2[:,:,0] - qq5 * dx4/kapp)
          w2[:,:,-1] =   w2[:,:,-1] + 2 *kapp* (dt4/(dx4**2)) * (w2[:,:,-2] - w2[:,:,-1] + qq6 * dx4/kapp)
         
      w2[1:,:-1,:-1] = np.nan    #We'll only plot the &quot;outside&quot; points.
      w2_uno = np.reshape(w2,-1)    
  
      #Plotting
      fig = pyplot.figure()
      X4, Y4, Z4 = np.meshgrid(x4, y4,z4)
      ax = fig.add_subplot(111, projection='3d')
      img = ax.scatter(X4, Y4, Z4, c=w2_uno, cmap=pyplot.jet())
      fig.colorbar(img)
      pyplot.show()
  
  </code></pre>
  <p>For 5000 iterations (qq6 = 1000/Area) we get:</p>
  <p><img src="https://i.imgur.com/3sBcgua.png" alt="" /></p>
  <p>I'm only applying heat on the top surface. Somehow I end up with the bottom surface heating up.</p>
  <p>Adding to that I'm trying to confine the region to which I apply heat (just a small part of one face). When I try that it seems like heat is only transferred in one direction, ignoring all others. Applying (qq1 = 1000/Area) for half of the front face (the other part is adiabatic, q = 0) we end up with:</p>
  <p><img src="https://i.imgur.com/J3TOkYm.png" alt="" /></p>
  <p>This is pretty odd. I suspect I'm having some trouble in the main loop (or maybe in the boundary conditions) that I'm not finding.</p>