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 [3]:
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[3]:
[<matplotlib.lines.Line2D at 0x7fc89533d5a0>]

Get some derivatives of the curve:

In [4]:
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[4]:
[<matplotlib.lines.Line2D at 0x7fc87f19e1d0>]

Get normals to the curve:

In [5]:
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[5]:
<matplotlib.quiver.Quiver at 0x7fc87f1c7af0>
In [6]:
dist_vec = tgt[:, np.newaxis] - path

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

pt.plot(dist)
Out[6]:
[<matplotlib.lines.Line2D at 0x7fc87ef93b80>]

Single-layer potential¶

Let's look at the integrand for the SLP:

In [7]:
slp_integrand = np.log(dist)
pt.plot(slp_integrand)

/tmp/ipykernel_384510/797241706.py:1: RuntimeWarning: divide by zero encountered in log
  slp_integrand = np.log(dist)
Out[7]:
[<matplotlib.lines.Line2D at 0x7fc87effe800>]

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

In [8]:
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])

/tmp/ipykernel_384510/4173338573.py:3: RuntimeWarning: divide by zero encountered in log
  log_sin_squared = 0.5*np.log(4*np.sin((tgt_t - t)/2)**2)
Out[8]:
[<matplotlib.lines.Line2D at 0x7fc87ee6d8a0>]
In [9]:
slp_subtracted = slp_integrand - log_sin_squared
pt.plot(slp_subtracted[near_sing_slice])

/tmp/ipykernel_384510/3811592322.py:1: RuntimeWarning: invalid value encountered in subtract
  slp_subtracted = slp_integrand - log_sin_squared
Out[9]:
[<matplotlib.lines.Line2D at 0x7fc87eed8c40>]
In [10]:
pt.plot(slp_subtracted)
Out[10]:
[<matplotlib.lines.Line2D at 0x7fc87ed426e0>]

How does this help?

Double-layer potential¶

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

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

/tmp/ipykernel_384510/4139646203.py:1: RuntimeWarning: invalid value encountered in true_divide
  grad_slp = dist_vec/dist**2
Out[11]:
[<matplotlib.lines.Line2D at 0x7fc87eda9990>]

S'¶

In [12]:
sp_integrand = np.sum(grad_slp * normals[:, tgt_index, np.newaxis], axis=0)
pt.plot(sp_integrand)
Out[12]:
[<matplotlib.lines.Line2D at 0x7fc87ec20b80>]

Questions¶

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