MathOverflow will be down for maintenance for approximately 3 hours, starting Monday evening (06/24/2013) at approximately 9:00 PM Eastern time (UTC-4).

## Finding the length of a cubic B-spline

Can I find an analytical solution to the the length of an 2-dimensional cubic B-spline? All I can find are chorded approximations and the opinion that the analytic solution is "unbearably gruesome". However, I believe if i had the solution to this elliptic integral: $\int_0^u \sqrt{\left(b_{1} + 2 b_{2} t + 3 b_{3} t^{2}\right)^{2} + \left(a_{1} + 2 a_{2} t + 3 a_{3} t^{2}\right)^{2}} dt$, I'd have my answer (at least on a per-segment basis). My goal is implement it in some software, so a little up-front grue is ok. Is it really not worth pursuing? I can't believe I'm better off with a numerical solution. I'm after 64-bit accuracy if that helps with the discussion.

Sympy cannot integrate it. I have not tried Math-CAD, Maxima, or other CASs. I do not have any significant math background (if that's not already obvious).

Edit (background info for Greg):

This spline will be used as a control-surface for conserved quantities in a physics analysis. Flux of these quantities must also be known. The spline must be second-derivative continuous and must pass through specified and an arbitrary number of knot points. A cubic B-spline seems right for the job. It is also required that the length of the spline or any arbitrary interval along the spline also be determinable to a near-machine-precision value. Unfortunately, this needs to be done quickly because the knots of the spline and the parametrization intervals are themselves being converged-upon to meet criteria which is dependent on these lengths. And to make matters worse, this is a time-transient solution (which means another level of iteration).

I am using Python to prototype this on a 64-bit multi-core PC. (I wince as I type 'Python' because it has a reputation for being slow.) However, I am leveraging numpy (for basic array math), scipy (for FITPACK, QUADPACK, etc.), and possibly sympy (for elliptical functions) for the heavy lifting. I plan to (if I can ever find an exact closed-form "gruesome" solution) to write it up in C or FORTRAN and compile it as it's own module.

It has to run on a PC as it is not a one-time-use algorithm. It will be used by others many times over in a design environment. (I guess that's another level of iteration again.)

-
Have you tried mathematica? Mathematica has all of the integrals. – Harry Gindi Feb 2 2010 at 9:40
Not all of the integrals because most functions cannot be integrated in closed form. – lhf Feb 2 2010 at 12:02
Here, all means more than any of the other CAS's. – Harry Gindi Feb 2 2010 at 12:30
If you know the roots of the polynomial inside the square root then Mathematica can give you a closed form in terms of elliptic functions. But I don't think this solution is of any use to anyone. Computing roots of polynomials is probably numerically trickier than computing integrals, and Mathematica gives a rather large expression. – Dan Piponi Feb 2 2010 at 17:42
Also, a flux calculation does not directly use the length of a curve. It is a different integral that can be expressed with a length factor in the integrand, but it is not the same numerical analysis problem as computing the length. It could be easier to compute flux accurately than length accurately. Machine precision is also a strange standard for splines, because the splines themselves typically approximate some other curve far short of machine precision. Your problem feels ill-posed and proper advice is at best hit-and-miss. – Greg Kuperberg Feb 4 2010 at 3:32

You should first accept the fact that it's an elliptic integral, and therefore doesn't have an elementary expression without elliptic functions. If you had a numerical library with elliptic functions, then great. Otherwise, you need to either implement elliptic functions yourself, or implement numerical integration of your integral.

I recommend numerical integration, just because in context it is conceptually simple and reliable. Your integrand has a fairly tame form: It can't blow up, the integrand is continuous, and the integrand is also real analytic unless it touches zero. In this situation, Gaussian integration has excellent properties. I don't feel like doing a precise calculation, but I would expect that for any choice of the coefficients, Gaussian quadrature with just 5 evaluation points already has to be fairly close to the exact answer for any choices of the coefficients.

The above is part of an answer, but not a complete answer you really want 64 bits of accuracy. Assuming that the integrand is real analytic, Gaussian quadrature or Clenshaw-Curtis will converge exponentially. It seems reasonable enough to use Clenshaw-Curtis, which lets you recycle evaluation points and has a predictable formula for the numerical integration weights, with more and more points until the answer looks accurate.

The only problem is in the annoying case in which the integrand touches zero, or comes close to touching zero, which can be interpreted geometrically as a point on the spline with nearly zero velocity. (Typically it looks like a cusp.) Then the integrand is NOT real analytic and these numerical methods do not converge exponentially. Or, in the near-zero case, the integrand is analytic but the exponential rate of convergence is slow. I'm sure that there are tricks available that will handle this case properly: You could cut out an interval near the bad point and do something different, or you could subtract off a known integral to tame the bad point. But at the moment I do not have specific advice for an algorithm that is both reasonable fast and reasonably convenient. Clenshaw-Curtis is convenient and usually very fast for this problem, but not all that fast in bad cases if you push it to 64 bits of precision.

Also, these methods can be thought of as a more sophisticated version of chordal approximation. Chordal approximation faces the same issue, except worse: It never converges at an exponential rate for a non-trivial cubic spline. If you want 64 bits, you might need a million chords.

Meanwhile, the GNU Scientific Library does have elliptic function routines. If you have elliptic functions, then again, your integral is not all that simple, but it is elementary. I don't know whether GSL or equivalent is available for your software problem. If it is, then an elliptic function formula is probably by far the fastest (for the computer, not necessarily for you).

In a recent comment, bpowah says "All I wanted to know is whether or not it was faster to compute the given integral numerically or exactly." Here is a discussion. Computing an integral, or any transcendental quantity, "exactly" is an illusion. Transcendental functions are themselves computed by approximate numerical procedures of various kinds: Newton's method, power series, arithmetic-geometric means, etc. There is an art to coding these functions properly. A competitive implementation of a function even as simple as sin(x) is already non-trivial.

Even so, I'm sure that it's faster in principle to evaluate the integral in question in closed form using elliptic functions. It could be hard work to do this right, because the first step is to factor the quartic polynomial under the square root. That already requires either the quartic formula (unfortunately not listed in the GNU Scientific Library even though it has the cubic) or a general polynomial solver (which is in GSL but has unclear performance and reliability). The solution also requires elliptic functions with complex arguments, even though the answer is real. It could require careful handling of branch cuts of the elliptic functions, which are multivalued. With all of these caveats, it doesn't seem worth it to work out an explicit formula. The main fact is that there is one, if you have elliptic functions available but not otherwise.

The merit of a numerical integration algorithm such as Gaussian quadrature (or Clenshaw-Curtis, Gauss-Kronrod, etc.) is that it is vastly simpler to code. It won't be as fast, but it should be quite fast if it is coded properly. The only problem is that the integrand becomes singular if it reaches 0, and nearly singular if it is near 0. This makes convergence much slower, although still not as slow as approximation with chords. With special handling of the near-singular points, it should still be fine for high-performance numerical computation. For instance, a polished strategy for numerical integration might well be faster than a clumsy evaluation of the relevant elliptic functions.

-
I was under the impression that he's just trying to compute something, not implement computation in a program he's writing. – Harry Gindi Feb 2 2010 at 9:40
As I read his question, he tried to symbolically integrate the expression (even though it's an elliptic integral) in the hope of using the formula in his code. – Greg Kuperberg Feb 2 2010 at 15:01
Thank you for for the insight and the well thought out response(s). I implemented the Clenshaw-Curtis quadratic approximation, but it was quite slow for my needs. I do have the elliptic function routines (libraries) at my disposal, within the software, however, as this is my first exposure to elliptic functions, I am at a loss. I'm eager to learn and determined to find a closed-form solution. I suspect it is naïve of me to think I could learn to do this on my own in a matter of days. If not, do you recommend a resource or an example that is similar? – bpowah Feb 3 2010 at 0:33
Clenshaw-Curtis is not a "quadratic" approximation, and I'm a little surprised that you found it "slow". You should say a little more about your speed and precision requirements for me to really understand the question. (Possibly including language and typical hardware.) I agree though that elliptic functions could be the best way to do this --- once I understand the real question. – Greg Kuperberg Feb 3 2010 at 6:53
Whoops. I meant to say quadrature, not quadratic. I'm using QUADPACK, but with the high precision I need and the fact that i need to generate many, many splines sequentially (inside another minimization routine), It does appear to be too slow for my needs. I'll post more info in my question later today. – bpowah Feb 4 2010 at 0:51

For numerical integration methods, see http://groups.google.com/group/comp.graphics.algorithms/msg/c7025fd53b18db94.

One of the links there is broken. Try "Adaptive sampling of parametric curves".

-

I presume you found this already and were not satisfied? As has been mentioned by other posters, you will still have to do things numerically, whether you use a quadrature rule or use any of the special algorithms for computing elliptic integrals.

-