Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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: enter image description here

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?

share|improve this question
    
@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
add comment

1 Answer 1

up vote 1 down vote accepted

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.

Maple craziness

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

share|improve this answer
    
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
show 2 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.