Kussmaul-Martensen quadrature (also often called "Kress quadrature")

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'?