Hacker News Reader

Comment by xyzzyz | Hacker News Reader

Why don’t you try the analytic method for cubics? The Abel-Ruffini only says that the anti-derivative is not elementary, not that it doesn’t exist. The eliptic integrals are non-elementary, but calculating them on a computer is not really different in principle than e.g. computing natural logarithm.

If you plug the desired integral into Mathematica or Wolfram Alpha, it should just spit out the relevant anti derivative expressed in terms of some standard elliptic function, which can then be calculated using some series expansion (or some library of math functions).

View on HN

You're really trying to nerd snipe, aren't you? I spent some quality time with elliptic functions as the solutions to the elastica, and also with the Fresnel integral which of course generates the Euler spiral. So yes, it wouldn't be surprising if there were a closed form solution, and that it might even be more efficient to compute.

So, this would be a good topic for someone else to show off their math skills :)

Haha, of course it wouldn’t really be much of an improvement over what you already have, and clearly I’m too lazy to do it myself — I would just enjoy looking at it :)

For those unfamiliar with the term 'nerd sniping' (as was I): https://xkcd.com/356/

Given that the desired integral will be a sixth order polynomial, it's not even that "we can't find the intergral" (we can), it's that there is no closed form solution once we've found it.

  Integrate[Sqrt[a + b u + c u^2 + d u^3], u]

Produces a huge expression which is a sum of EllipticF[ArcSin[#]] terms. I think it would be better to try and numerically compute the integral.

No, that’s precisely my point: this expression can be evaluated faster and more accurately than by applying generic numeric integration method. Computing this EllipticF function is not different than computing natural logarithm, and if you don’t mind latter, there is not much reason to shy from former.

To be specific

r1,r2,r3 are roots of a+bu+cu^2+u^3,

Cleaning up the output a bit, you get

(2 (b c EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)]-9 a d EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)]+(c+3 d u) (a+u (b+u (c+d u))) Sqrt[((-u+r2) (u-r3))/(r2-r3)^2] (r2-r3)+2 c^2 (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)] (EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] r1+EllipticE[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (-r1+r3))-6 b d (u-r2) (u-r3) Sqrt[(-u+r1)/(r1-r3)] (EllipticF[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] r1+EllipticE[ArcSin[Sqrt[(-u+r3)/(-r2+r3)]],(r2-r3)/(r1-r3)] (-r1+r3))))/(15 d Sqrt[a+u (b+u (c+d u))] Sqrt[-(((u-r2) (u-r3))/(r2-r3)^2)] (r2-r3)).

But I guess if EllipticF is fast to compute, finding the roots can also be done pretty quickly, so yes its probably not that bad.

You need to integrate the square root of a quartic polynomial – the derivative of a cubic is a quadratic, then you square the norm of that, which gives a quartic. I spent a little time in Mathematica and it there is a solution, but it seems very unlikely to me it's going to be more efficient than numerical integration - it involves finding the roots of a quartic and has nine elliptic function evaluations. From a slide I found[0] it looks like it takes on the order of 100ns to compute an elliptic function, so at best it's still going to be slower than the adaptive quadrature even with an accuracy of 1e-12 or so. Plus there are a lot of divisions by differences of roots, and as each one of those approaches zero there will be numerical stability issues.

Long story short, it looks doable, but is a lot of work to get right and has almost no chance of beating out numerical integration in performance.

[0] http://arith22.gforge.inria.fr/slides/05-fukushima.pdf

Reply on news.ycombinator.com