Inaccuracy in numerical calculation of arclength of part of an ellipse

Question

I am trying to numerically calculate the arclength of part of an ellipse according to:

$$ L = \int_0^{\phi_s}\sqrt{r^2+\left(\frac{dr}{d\phi}\right)^2} d\phi $$

where $r$ is defined as: $$ r=\frac{a b}{\sqrt{a^2\sin^2\phi+b^2\cos^2\phi}} $$ with $2a$ and $2b$ the major and minor axis of the ellipse

This seems to work well in general, but it goes wrong if $\phi_s$ is small and $a$ and $b$ are big. In that case the calculation is horribly off ($\approx$ factor 5).

The specific example I am working on has $\phi_s=2.06^\circ (=0.036 \text{ rad})$, $a=16.3$ and $b=0.52$ and I find $L=5.48$. Judging by distance between the points $\phi_s=0$ and $\phi_s=2.06$ which are located at relative distances (0.572,0.430) on the ellipse I would expect a value for $L$ close to 0.8 instead.

Interestingly, I checked and the numerical integration itself is not the issue! If I plot the integrant and roughly compute $L$ I would indeed find 5.5 as you can see in the plot below:

Can anybody explain what the problem is? Is this related to round-off errors somehow? Or do I get some form of catastrophic cancellation somewhere?! Because I have the problem both in mathematica and in matlab I wonder whether I am running into machine-precision issues somehow?

Answer 1

This is what I get from Maple 16.

a:=16.3; b:=0.52; phis:=0.036; 
r:=a*b/sqrt(a^2*sin(phi)^2+b^2*cos(phi)^2);
L:=Int(sqrt(r^2+diff(r,phi)^2),phi=0..phis);
evalf(L)

Answer: 5.526897611 which I think is correct.

According to Maple, the integrand is $$ab \sqrt{\frac{b^4\cos^2\phi+a^4\sin^2\phi}{(b^2\cos^2\phi+a^2\sin^2\phi)^3}}$$ (found using simplify(sqrt(r^2+diff(r,phi)^2)); without entering numerical values for $a,b$, and then simplifying further by hand).

---

Added. This may be unrelated, but is a curious observation anyway. Maple 16 creates a messy formula if I ask it to simplify sqrt(r^2+diff(r,phi)^2) after giving the numerical values, such as a:=16.3; b:=2.5; Here is the screenshot.

Unsurprisingly, an attempt to use this "simplified" expression chokes up both plotting and numerical integration routines.