In [1]:
import numpy as np
import matplotlib.pyplot as pt

In [2]:
t = np.linspace(0, 2*np.pi, 300,endpoint=False)


## Setup¶

Let's make a curve and pick a target point:

In [29]:
uncircleness = 1

path = np.array([
np.cos(t) + uncircleness*0.2*np.sin(3*t),
np.sin(t) + uncircleness*0.1*np.sin(3*t)
])

tgt_index = len(t)//2
tgt_t = t[tgt_index]
tgt = path[:, tgt_index]

pt.gca().set_aspect("equal")
pt.plot(path[0], path[1])
pt.plot(tgt[0], tgt[1], "o")

Out[29]:
[<matplotlib.lines.Line2D at 0x7f55180f0048>]

Get some derivatives of the curve:

In [28]:
import scipy.fftpack as fft

dpath_dt = np.array([
fft.diff(path[0]),
fft.diff(path[1]),
])

dpdt_squared = dpath_dt[0]**2 + dpath_dt[1]**2
pt.plot(dpdt_squared)

Out[28]:
[<matplotlib.lines.Line2D at 0x7f551810a668>]

Get normals to the curve:

In [27]:
normals = np.array([
dpath_dt[1],
-dpath_dt[0]
]) / np.sqrt(dpdt_squared)

pt.plot(path[0], path[1])
pt.quiver(path[0], path[1], normals[0], normals[1])

Out[27]:
<matplotlib.quiver.Quiver at 0x7f55181c4c18>
In [26]:
dist_vec = tgt[:, np.newaxis] - path

dist = np.sqrt(np.sum(dist_vec**2, axis=0))

pt.plot(dist)

Out[26]:
[<matplotlib.lines.Line2D at 0x7f55181b4ef0>]

## Single-layer potential¶

Let's look at the integrand for the SLP:

In [25]:
slp_integrand = np.log(1/dist)
pt.plot(slp_integrand)

-c:2: RuntimeWarning: divide by zero encountered in true_divide

Out[25]:
[<matplotlib.lines.Line2D at 0x7f551824f6d8>]

Even if this is integrable--Gaussian quadrature will do a terrible job. Why?

In [24]:
near_sing_slice = slice(tgt_index-20, tgt_index+20)

log_sin_squared = -0.5*np.log(4*np.sin((tgt_t - t)/2)**2)
pt.plot(log_sin_squared[near_sing_slice])
pt.plot(slp_integrand[near_sing_slice])

-c:4: RuntimeWarning: divide by zero encountered in log

Out[24]:
[<matplotlib.lines.Line2D at 0x7f5518315b70>]
In [23]:
slp_subtracted = slp_integrand - log_sin_squared
pt.plot(slp_subtracted[near_sing_slice])

-c:2: RuntimeWarning: invalid value encountered in subtract

Out[23]:
[<matplotlib.lines.Line2D at 0x7f55182f6358>]
In [17]:
pt.plot(slp_subtracted)

Out[17]:
[<matplotlib.lines.Line2D at 0x7f55184f2e80>]

How does this help?

### Double-layer potential¶

In [22]:
grad_slp = dist_vec/dist**2

dlp_integrand = np.sum(grad_slp * normals, axis=0)
pt.plot(dlp_integrand)

-c:2: RuntimeWarning: invalid value encountered in true_divide

Out[22]:
[<matplotlib.lines.Line2D at 0x7f551838fb38>]

### S'¶

In [21]:
sp_integrand = np.sum(grad_slp * normals[:, tgt_index, np.newaxis], axis=0)
pt.plot(sp_integrand)

Out[21]:
[<matplotlib.lines.Line2D at 0x7f55183b13c8>]

## Questions¶

• How would you apply this for Helmholtz?
• Name aspects that make this rule slightly impractical
• How would this apply to D'?