1. import numpy as np
  2.  
  3. nx=3
  4. nz=4
  5. flux2 = np.empty((nx+1, nz+1))
  6. flux2[:] = np.random.normal(size=(nx+1, nz+1))
  7. flux3 = flux2.copy()
  8. bz = np.random.normal(size=(nx+1, nz+1))
  9. bx = np.random.normal(size=(nx+1, nz+1))
  10.  
  11. dz = 0.112344
  12. dx = 0.015543
  13. for i in np.arange(1,nx+1):
  14.     for j in np.arange(1,nz+1):
  15.         flux2[i,j] = 0.5*(flux2[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux2[i,j-1] - bx[i,j-1]*dz)
  16.  
  17. a = 0.5
  18. aexp = np.arange(nz).reshape(nz, 1) - np.arange(nz).reshape(1, nz)
  19. abcoeff = 0.5**aexp
  20. abcoeff[aexp<0] = 0
  21. for i in np.arange(1,nx+1):
  22.     b = 0.5*flux3[i-1, 1:] + 0.5*bz[i-1, 1:]*dx - 0.5*bx[i,:-1]*dz
  23.     bvals = (abcoeff * b.reshape(1, nz)).sum(axis=1)
  24.     n = np.arange(1, nz+1)
  25.     x0 = flux3[i, 0]
  26.     #ab = b * np.cumsum(a**(n-1))
  27.     flux3[i, 1:] = a**n * x0 + bvals
  28.  
  29. print np.abs(flux2-flux3) < 1e-15