Incomplete Integrals

[zwhitchcox]

Could you please add the incomplete elliptic integral of the second kind? Or is that too difficult? Trying to figure it out now and can't figure out how to do it.

[duetosymmetry]

Hi Zane, it's a good idea, but I haven't spent much time digging into the literature and thinking about implementing this.

My understanding is that the modern approach is to convert classic elliptic integrals into Carlson symmetric form, and just implement the Carlson integrals.

Please let me know if you find pseudocode listings for how to implement the Carlson symmetric integrals robustly.

[duetosymmetry]

The efficient way to implement these is discussed in Carlson's proceedings Elliptic Integrals: Symmetry and symbolic integration [PDF], from the 1997 Turin conference on the centennial of Tricomi's birth. The duplication formulas are summarized at e.g. https://dlmf.nist.gov/19.26#iii. See also Carlson's arXiv:math/9409227 and Johansson's arXiv:1806.06725.

[tibotiber]

Hi, I landed here as I needed incomplete integrals to compute equidistant points on an ellipse (same arc length). In case this can help anyone else, I came across this numerical integration method which is not exact but was good enough for my use-case: https://stackoverflow.com/a/6972434. Here is my JS implementation:

export const getEquidistantPointsOnEllipse = ({
  r1,
  r2,
  numberOfPoints,
  deltaTheta = degToRad(1)
}) => {
  // we want to avoid calculating the inverse of the incomplete elliptic integral of the second kind
  // https://math.stackexchange.com/questions/172766/calculating-equidistant-points-around-an-ellipse-arc
  // could only find complete elliptic integrals in JS
  // https://github.com/duetosymmetry/elliptic-integrals-js/issues/2
  // we trick by integrating the arc perimeter manually and marking the points as we do
  // https://stackoverflow.com/questions/6972331/how-can-i-generate-a-set-of-points-evenly-distributed-along-the-perimeter-of-an
  // # dp(t) = sqrt( (r1*sin(t))^2 + (r2*cos(t))^2)
  // # circ = sum(dp(t), t=0..2*Pi step 0.0001)
  // # n = 20
  // # nextPoint = 0
  // # run = 0.0
  // # for t=0..2*Pi step 0.0001
  // #     if n*run/circ >= nextPoint then
  // #         set point (r1*cos(t), r2*sin(t))
  // #         nextPoint = nextPoint + 1
  // #     next
  // #     run = run + dp(t)
  // # next
  const integralSteps = Math.round((Math.PI * 2) / deltaTheta)
  const dp = t =>
    Math.sqrt(Math.pow(r1 * Math.sin(t), 2) + Math.pow(r2 * Math.cos(t), 2))
  const circ = compose(
    sum,
    map(step => dp(step * deltaTheta))
  )(range(0, integralSteps))
  // then gather equi-arcdistant points along the ellipse
  let run = 0
  let nextPoint = 0
  let points = []
  for (let i = 0; i < integralSteps; i++) {
    const t = i * deltaTheta
    if ((numberOfPoints * run) / circ >= nextPoint) {
      points.push({
        radAngle: t,
        x: r1 * Math.cos(t),
        y: r2 * Math.sin(t)
      })
      nextPoint++
    }
    run += dp(t)
  }
  return points
}

[duetosymmetry]

Hi @tibotiber, thanks for the message. I am open to implementing additional functions like incomplete integrals, their inverses (for example the Jacobi sn() function), and related functions (like the Jacobi am() function). I would avoid algorithms that rely on direct quadrature because they do not have the convergence guarantees that iterative ones enjoy (e.g. based on duplication formulas or the AGM), rather they converge with the step size of the quadrature routine, making their speed depend on the argument.

I would suggest that you take a look at https://en.wikipedia.org/wiki/Jacobi_elliptic_functions#Definition_as_trigonometry:_the_Jacobi_ellipse and see if you can rephrase your problem in terms of am(), which enjoys a rapidly converging iterative algorithm, which I've seen in the literature.

[forresto]

I ended up here b/c I'm trying to translate this drawing to JS...
https://loopspace.mathforge.org/CountingOnMyFingers/Triskele/#section.4

They have python code that includes from scipy.special import ellipe, ellipeinc

e = np.sqrt(1 - args.height**2)
a = .5

for i in range(args.steps+1):
    phi = float(i)/args.steps*np.pi - np.pi/2
    fd.write("{} {}\n".format(args.height * a * np.cos(phi), a * (ellipeinc(phi,e) + ellipe(e))))