title: Apply dynamical function to each point in phase space (represented by a 2D matrix)
body: <p>I have a matrix of integers, <code>phase_space</code> of shape <code>(n,n)</code>, where each entry represents the number of points in that location in space. I also have two update matrices <code>u_x, u_y</code> also of shape <code>(n,n)</code>, with integers in the range <code>0,...,n</code> specifying where my dynamical system takes each corresponding point in space.
  I want to "apply" the update matrices to the phase space iteratively.</p>
  
  <p>For example, if</p>
  
  <pre><code>&gt;&gt;&gt;u_x
  array([[1, 2, 1],
         [0, 1, 2],
         [0, 0, 0]])
  &gt;&gt;&gt;u_y
  array([[2, 1, 2],
         [1, 0, 1],
         [2, 2, 0]])
  &gt;&gt;&gt;phase_space 
  array([[1, 1, 1],
         [1, 1, 1],
         [1, 1, 1]])
  </code></pre>
  
  <p>I want</p>
  
  <pre><code>&gt;&gt;&gt;new_phase_space
  array([[1., 1., 2.],
         [1., 0., 2.],
         [0., 2., 0.]])
  </code></pre>
  
  <p>My current (working) solution is to loop as follows</p>
  
  <pre><code>for i in range(n):
      for j in range(n):
          new_phase_space[u_x[i, j], u_y[i, j]] += phase_space[i,j] 
  </code></pre>
  
  <p>Is there any way to vectorize this?</p>