# 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)

Out[32]:
<matplotlib.image.AxesImage at 0x7fb1b84b5cc0>


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()

Out[75]:
<matplotlib.colorbar.Colorbar at 0x7fb1b06590b8>

• 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()

Out[99]:
<matplotlib.colorbar.Colorbar at 0x7fb1a99dbc18>

• 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)

Out[115]:
<matplotlib.image.AxesImage at 0x7fb1a95546a0>


Why is this thing called a 'multipole' expansion?