Saturday, January 15, 2011

Discrete Laplace Transform - del2

I ported Matlab del2() call into Python. Here it is for reference:
`import numpy as npdef del2(M):  dx = 1  dy = 1  rows, cols = M.shape  dx = dx * np.ones ((1, cols - 1))  dy = dy * np.ones ((rows-1, 1))  mr, mc = M.shape  D = np.zeros ((mr, mc))  if (mr >= 3):      ## x direction      ## left and right boundary      D[:, 0] = (M[:, 0] - 2 * M[:, 1] + M[:, 2]) / (dx[:,0] * dx[:,1])      D[:, mc-1] = (M[:, mc - 3] - 2 * M[:, mc - 2] + M[:, mc-1]) \          / (dx[:,mc - 3] * dx[:,mc - 2])         ## interior points      tmp1 = D[:, 1:mc - 1]      tmp2 = (M[:, 2:mc] - 2 * M[:, 1:mc - 1] + M[:, 0:mc - 2])      tmp3 = np.kron (dx[:,0:mc -2] * dx[:,1:mc - 1], np.ones ((mr, 1)))      D[:, 1:mc - 1] = tmp1 + tmp2 / tmp3  if (mr >= 3):      ## y direction      ## top and bottom boundary      D[0, :] = D[0,:]  + \          (M[0, :] - 2 * M[1, :] + M[2, :] ) / (dy[0,:] * dy[1,:])         D[mr-1, :] = D[mr-1, :] \          + (M[mr-3,:] - 2 * M[mr-2, :] + M[mr-1, :]) \          / (dy[mr-3,:] * dx[:,mr-2])         ## interior points      tmp1 = D[1:mr-1, :]      tmp2 = (M[2:mr, :] - 2 * M[1:mr - 1, :] + M[0:mr-2, :])      tmp3 = np.kron (dy[0:mr-2,:] * dy[1:mr-1,:], np.ones ((1, mc)))      D[1:mr-1, :] = tmp1 + tmp2 / tmp3     return D / 4`