# Inaccuracy in numerical calculation of arclength of part of an ellipse

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?

-
@your first comment: sorry I wrote down the value in degrees instead of radian, I fixed it in the post. The relative distance comes from the fact that I know where the tip of the ellipse and the point $\phi_s$ are ---- your second comment: So you are saying that the problem is most likely caused by the parametrization in combination with round-off errors? – Michiel May 3 '13 at 6:00
yeah, I see that now. Do you have a symbolic answer to the integrand? Then at least I can start looking for my mistake from there! Anyways, Thanks for the quick help!! – Michiel May 3 '13 at 6:03
By the way: I agree that a different parametrization would be better. The thing is that I also need the polar angle so I will have to do a bit of work connecting the parametrization you suggest with the polar angle – Michiel May 3 '13 at 6:06

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.

-
mmm that's odd. I found this for the integrand: $$a b \sqrt{\frac{b^4 \cos^2\phi+a^4 \sin^2\phi}{(b^2\cos^2\phi+a^2\sin^2\phi)^3}}$$ which is supposed to be the same according to ($\cos^2\phi+\sin^2\phi=1$) but apparently it is not???? – Michiel May 3 '13 at 6:14
@michielm Actually it's the same thing: $1-\cos^2 = \sin^2$; Maple (and myself) were too lazy to simplify. – 75064 May 3 '13 at 6:16
@michielm I can't believe it. Any sensible program will handle both forms in the same way. I get 0.74... by integrating the formula in your comment. – 75064 May 3 '13 at 6:22
You're right, they both are the same (I was missing a bracket), but still I get a much higher value. What I noticed by the way is that you use b=2.5 while I have b=0.52. What do you get with that smaller b? – Michiel May 3 '13 at 6:35
@michielm Hm, indeed I had a different $b$. With $b=0.52$ I get $5.52$ from the integral. But this looks right to me. The $r$-value drops from 16.3 at $\phi=0$ to about $10.8$ at $\phi=0.036$. So the arc could not be any shorter than $16.3-10.8$. The ellipse is practically a needle. – 75064 May 3 '13 at 6:35