# coding: utf-8 # # Multipole and Local Expansions # In[1]: import numpy as np import numpy.linalg as la import matplotlib.pyplot as pt # Let's consider a potential. This one could look slightly familiar from a homework assignment. # In[31]: nsources = 15 nxtgts = 400 nytgts = 400 angles = np.linspace(0, 2*np.pi, nsources, endpoint=False) r = 1 + 0.3 * np.sin(3*angles) sources = np.array([ r*np.cos(angles), r*np.sin(angles), ]) charges = 0.5+2*np.sin(3*angles) left, right, bottom, top = extent = (-2, 4, -4, 2) targets = np.mgrid[left:right:nxtgts*1j, bottom:top:nytgts*1j] # In[32]: pt.plot(sources[0], sources[1], "x") dist_vecs = sources.reshape(2, -1, 1, 1) - targets.reshape(2, 1, targets.shape[-1], -1) dists = np.sqrt(np.sum(dist_vecs**2, axis=0)) potentials = np.sum(charges.reshape(-1, 1, 1) * np.log(dists), axis=0) pt.imshow(potentials.T[::-1], extent=extent) # Now let's create a stash of derivatives, all about a center of 0, to make things easier: # In[33]: def f(arg): return np.log(np.sqrt(np.sum(arg**2, axis=0))) def fdx(arg): x, y = arg r2 = np.sum(arg**2, axis=0) return x/r2 def fdy(arg): x, y = arg r2 = np.sum(arg**2, axis=0) return y/r2 def fdxx(arg): x, y = arg r2 = np.sum(arg**2, axis=0) return 1/r2 - 2*x**2/r2**2 def fdyy(arg): x, y = arg r2 = np.sum(arg**2, axis=0) return 1/r2 - 2*y**2/r2**2 def fdxy(arg): x, y = arg r2 = np.sum(arg**2, axis=0) return - 2*x*y/r2**2 # ## Local expansions # In[71]: center = np.array([1.5, -1]) #center = np.array([2, -2]) #center = np.array([3, -3]) #center = np.array([0, 0]) # In[72]: expn = 0 for isrc in range(nsources): a = sources[:, isrc] - center hx, hy = center.reshape(-1, 1, 1) - targets expn += charges[isrc]*( f(a) + fdx(a)*hx + fdy(a)*hy + fdxx(a)*hx**2/2 + fdxy(a)*hx*hy + fdxy(a)*hy**2/2 ) # In[75]: err = expn - potentials pt.plot(center[0], center[1], "o") pt.plot(sources[0], sources[1], "x") pt.imshow(np.log10(1e-2 + np.abs(err.T[::-1])), extent=extent) pt.colorbar() # * Move the center around, see how the errors change # * Reduce to linears, see how the errors change # ## Multipole expansions # In[90]: center = np.array([0, 0]) # center = np.array([1, 0]) # In[101]: expn = 0 for isrc in range(nsources): a = targets - center.reshape(-1, 1, 1) hx, hy = center - sources[:, isrc] expn += charges[isrc]*( f(a) + fdx(a)*hx + fdy(a)*hy + fdxx(a)*hx**2/2 + fdxy(a)*hx*hy + fdxy(a)*hy**2/2 ) # In[99]: err = expn - potentials pt.plot(center[0], center[1], "o") pt.plot(sources[0], sources[1], "x") pt.imshow(np.log10(1e-2 + np.abs(err.T[::-1])), extent=extent) pt.colorbar() # * Move the center around, observe convergence behavior # * Reduce to linears, observe convergence behavior # Look at individual basis functions: # In[115]: pt.imshow( fdx(targets).T[::-1], extent=extent, vmin=-1, vmax=1) # Why is this thing called a 'multipole' expansion?