I'm using a cubic spline interpolation for given data points. The boundary condition for the spline is that $f'(a)$ and $f'(b)$ are given (I'm using a finite difference formula $\frac{y_1-y_0}{x_1-x_0}$ to calculate the first derivatives $f'(a)$ and $f'(b)$, so they are not really given). With $a$ and $b$ being the boundaries. I ran into the problem, that my second derivative becomes really inaccurate at the boundary, if $f'(a)$ and $f'(b)$ aren't good approximations. I can't use a natural spline, since I also need a value for the second derivative at the endpoints. I'm using the spline to solve a stiff PDE. Are there any other good methods to approximate the first derivative at the endpoints?
Good Approximation for First Derivative at Endpoints for Cubic Spline Interpolation
differential equationsna.numerical-analysissplines
Related Solutions
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.
If there are no more constraints, then you can do it with arbitrarily low curvature with any reasonable class of splines. If the points are say within a 10cm region, make huge loops 1km in diameter (or bigger if you want smaller curvature). If the spline construction is smooth, continuous, and invariant under similarity, then the curvature converges to 0.
If the curves are required to stay in a bounded region of the plane, then as the region gets smaller, not even arbitrary $C^2$ curves can thread through them with bounded curvature. Just imagine $3$ points at the corners of an equilateral triangle 1 micron on a side, and ask for the curve to be confined to a box 2 microns on a side. The curvature will be on the order of $\pi / $ micron.
Here are two copies of a set of four points threaded with Adobe Illustrator splines to illustrate the phenomenon. Note: I added extra knots in the big loops to make them look better, but this isn't necessary to construct examples. (The mathematical characterization of these splines is not relevant to the answer, and furthermore, I don't actually know):
alt text http://dl.dropbox.com/u/5390048/splines.jpg
The design considerations for splines are much more subtle than minimizing curvature.
However, I'd like to mention that the earlier meaning of splines had to do with thin splints of wood used in woodworking, e.g. boat building, to lay out curves for cutting. Unlike the usual mathematical splines, they have fixed length, and a reasonable mathematical model is that they trace out curves that locally minimize total curvature subject to their constraints (lead weights called ducks because that's what they resembled).
It's easy to get examples of these traditional splines with multiple local minima: cut a strip of paper (good enough for this) and bring the ends closer without turning them. The strip pops to one side or the other, giving two local minima.
Best Answer
For instance, you can approximate $f(x)$ by the Newton interpolating polynomial $P(x)$ based on a number of nodes closest to (and including) a given endpoint $x_0$, and then approximate $f'(x_0)$ by $P'(x_0)$.
Here is an image of a Mathematica notebook illustrating this approach on simulated data: