Problem: I am currently trying to integrate a singular kernel function of the type $$G(x,y)=\frac{\exp(ik||x-y||_2)}{4\pi ||x-y||_2}$$ which lies at the centre of a triangle, over this triangle. $i$ is the imaginary unit, $k$ denotes the wavenumber. There is some literature on this in the context of potential problems such as the "singularity subtraction method" but I can't find a simple code example in Matlab or Python.

Question: Is there anyone who implemented this or could point me to simple matlab/python code for a reference triangle or to a scipy routine or similar that performs this numerical integration?

  • $\begingroup$ What exactly is the integral you are trying to compute? What lies at the center of the triangle? Which of the 52,440 triangle centers? $\endgroup$ Commented Feb 2, 2023 at 16:28
  • 3
    $\begingroup$ Please see my answer for the EFIE (electric field integral equation) question. We can close this one as a duplicate, or you may want to edit it to point out a difference. $\endgroup$
    – Anton Menshov
    Commented Feb 2, 2023 at 17:29
  • $\begingroup$ Please do not close the question as I explicitly asked for a simple example code and it is not mentioned in the answers of this question or the question you linked. $\endgroup$
    – Bulbasaur
    Commented Feb 3, 2023 at 7:43

1 Answer 1


This sort of integral often arises in the context of surface integral equation methods for the Helmholtz or Maxwell PDE's. A good entry-level method to cancel the singularity is the "Duffy transform", which splits the triangle at the singularity into 3 subtriangles (each with a singularity at a vertex), then uses a clever change/transform of coordinate variables, to rewrite in the integral in such a way that the Jacobian of the transformation tidily cancels the singularity. There's actually a lot of methods in this space, searching for "singularity cancellation" or "singularity subtraction" will yield numerous results.

If I recall correctly, a useful trick in the context of some Galerkin/method-of-moments formulations, is to split the numerator $exp(..)$ into two terms over the same denominator $d$, something like $exp(..)/d = [exp(..)-1]/d + 1/d$. The first term is actually well behaved / smooth, as its numerator vanishes to zero at the same rate as the denominator. It can be integrated accurately using Gaussian quadrature as long as you use distinct rules for the test and source loops (so that you avoid direct/naive evaluation at x=y, which would yield 0/0). The second term remains singular but is now purely real-valued, which can be a nice simplification for debugging and also reduces the underlying arithmetic costs during Duffy transformation (or whatever else you might select).

  • $\begingroup$ Do you know of a publication that shows Matlab/Python code or a routine where you can see what is happening in the code? $\endgroup$
    – Bulbasaur
    Commented Feb 3, 2023 at 7:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.