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

Copyright (C) 2020 Andreas Kloeckner

MIT License Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

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>]
No description has been provided for this image

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>]
No description has been provided for this image

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>
No description has been provided for this image
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>]
No description has been provided for this image

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>]
No description has been provided for this image

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>]
No description has been provided for this image
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>]
No description has been provided for this image
In [10]:
pt.plot(slp_subtracted)
Out[10]:
[<matplotlib.lines.Line2D at 0x7fc87ed426e0>]
No description has been provided for this image

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>]
No description has been provided for this image

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>]
No description has been provided for this image

Questions¶

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