The Surface Area Of An Ellipsoid

A. Dieckmann, Universität Bonn, July 2003

    This short note shows a way to the formula for the surface area of an ellipsoid. The result is given for example in the formula gallery of the Mathematica Book (Ref. 1), but with a typo. There was no instruction accessible to me, how it is obtained, and I had to try for myself. It proved harder than I first thought, so I think it is worthwhile to put the derivation on the internet.
    To make a later integration valid, we have to name the semi-axes of the ellipsoid c > a > b > 0, so that b < c and   0<(k=c^( 2)/(c^2-b^2)*(a^2-b^2)/a^2) <1. The axes may be renamed at the end.

a = 12/10 ; b = 11/10 ; c = 47/10 ; (c^2 (a^2 - b^2))/(a^2 (c^2 - b^2))//N

0.168978

Now calculate the area numerically – for reference:

NIntegrate[Sin[θ] (a^2 b^2 Cos[θ]^2 + b^2c^2 Sin[θ]^2 Cos[φ]^2 + a^2 c^2 Sin[θ]^2Sin[φ]^2)^(1/2), {θ, 0, π}, {φ, 0, 2π}]

54.6901

Symmetry allows the reduction of the integration domain:

54.69012410

Integrate symbolically with respect to θ:

a=. ; b=. ; c=. ;

Integrate[Sin[θ] (a^2 b^2 Cos[θ]^2 + b^2c^2 Sin[θ]^2 Cos[φ]^2 + a^2 c^2 Sin[θ]^2Sin[φ]^2)^(1/2), {θ, 0, π/2}] * 8//FullSimplify

After numerical integration with respect to φ the area comes out correctly:

a = 1.2 ; b = 1.1 ; c = 4.7 ;

54.6901 + 0. 

Express ArcSinh by ArcTanh:

54.6901 + 0. 

a=. ; b=. ; c=. ;

And the integral can be reduced to a simpler form by substituting  c^( 2) (Cos[φ]^( 2)/a^( 2) + Sin[φ]^( 2)/b^( 2) ) = x:

Solve[x == c^2 (Cos[φ]^2/a^2 + Sin[φ]^2/b^2), φ]

Solve :: ifun : Inverse functions are being used by Solve, so some solutions may not be found.

The Jacobian dφ/dx is:

D[ArcCos[(a (c^2 - b^2 x)^(1/2))/((a^2 - b^2)^(1/2) c)], x]//FullSimplify

(a b^2)/(2 ((a - b) (a + b))^(1/2) c (b^2 (c^2 - a^2 x))/((-a^2 + b^2) c^2)^(1/2) (c^2 - b^2 x)^(1/2))

Check the new, 'simple' form :

a = 1.2 ; b = 1.1 ; c = 4.7 ;

2π a b - 2a^2 b^2NIntegrate[( ArcTanh[(1 - x)^(1/2)] x)/((1 - x)^(1/2) (c^2 - a^2 x)^(1/2) (-c^2 + b^2 x)^(1/2)), {x, c^2/a^2, c^2/b^2}, MinRecursion→10, MaxRecursion→33]

54.6901 + 0. 

Then integrate by parts with  u=ArcTanh[(1 - x)^(1/2)] , v'= x/((1 - x)^(1/2)(c^( 2) - a^( 2) x)^(1/2)(-c^( 2) + b^( 2) x)^(1/2)).

a=. ; b=. ; c=. ;

u ' = D[ArcTanh[(1 - x)^(1/2)], x]

-1/(2 (1 - x)^(1/2) x)

v = Integrate[( x)/((1 - x)^(1/2) (c^2 - a^2 x)^(1/2) (-c^2 + b^2 x)^(1/2)), x]

a = 1.2 ; b = 1.1 ; c = 4.7 ; uv[c^2/b^2] - uv[c^2/a^2]

-13.1977 + 0. 

uv(upper limit) - uv(lower limit):

(2ArcTanh[(1 - c^2/a^2)^(1/2)])/(a (b^2 - c^2)^(1/2)) ((1 - c^2/b^2) EllipticE[(c^2 (b^2 - a^2))/(a^2 (b^2 - c^2))] - EllipticK[(c^2 (b^2 - a^2))/(a^2 (b^2 - c^2))])

-13.1977 + 0. 

Surface area= 2π a b - 2a^( 2)b^( 2)(u v-∫u '  v):

54.6901 + 0. 

There are now two integrals left to attack, but luckily, both have already been done – about 130 years ago. The first is:

NIntegrate[1/(x (1 - x)^(1/2)) EllipticE[ArcSin[(-c^2/b^2 + x)/(c^2/a^2 - c^2/b^2)^(1/2)], (c^2 (b^2 - a^2))/(a^2 (b^2 - c^2))], {x, c^2/a^2, c^2/b^2}]

0. - 0.0348228 

a=. ; b=. ; c=. ;

Substitute ArcSin [((-c^( 2)/b^( 2) + x)/(c^( 2)/a^( 2) - c^( 2)/b^( 2)))^(1/2)] = z:

Solve[ArcSin[(-c^2/b^2 + x)/(c^2/a^2 - c^2/b^2)^(1/2)] == z, x]

{{x→c^2 (1/b^2 + Sin[z]^2/a^2 - Sin[z]^2/b^2)}}

D[c^2 (1/b^2 + Sin[z]^2/a^2 - Sin[z]^2/b^2), z]//FullSimplify

2 (1/a^2 - 1/b^2) c^2 Cos[z] Sin[z]

a = 12/10 ; b = 11/10 ; c = 47/10 ;

-0.0348228 

If the conditions cited at the beginning
• 0 < ( t = ArcCos[b/c] ) < π/2, or equivalently b < c, and
• 0 < k < 1
are all fulfilled, then this integral has a solution, found with Gradshteyn Ryzhik (6.123) (Ref.2) to be :

((b - a) b π)/(a (b^2 - c^2)^(1/2)) + 2ArcTanh[(1 - c^2/a^2)^(1/2)] EllipticE[(c^2 (b^2 - a^2))/(a^2 (b^2 - c^2))] + π/IEllipticE[ArcCos[b/c], (c^2 (b^2 - a^2) )/(a^2 (b^2 - c^2))]//N

0. - 0.0348228 

The second integral is very similar to the first :

NIntegrate[1/(x (1 - x)^(1/2)) EllipticF[ArcSin[(-c^2/b^2 + x)/(c^2/a^2 - c^2/b^2)^(1/2)], (c^2 (b^2 - a^2))/(a^2 (b^2 - c^2))], {x, c^2/a^2, c^2/b^2}]

0. - 0.0364505 

-0.0364505 

This is – with Gradshteyn Ryzhik (6.113-2) – :

2ArcTanh[ ( 1 - c^2/a^2)^(1/2)] EllipticK[(c^2 (b^2 - a^2))/(a^2 (b^2 - c^2))] + π/IEllipticF[ArcCos[b/c], (c^2 (b^2 - a^2))/(a^2 (b^2 - c^2))]//N

0. - 0.0364505 

Put them into the formula for the surface area above :

54.690124099820308

This is already the correct expression. Now rename the axes, that a > b > c :

a=. ; b=. ; c=. ;

aa = b ; bb = c ; cc = a ;

The final result is ( incidentally valid for any set of semi-axes – if a ≠ b ≠ c) :

a = 4.7 ; b = 1.2 ; c = 1.1 ; t = ArcCos[c/a] ; s = ArcCos[c/b] ; k = Sin[s]/Sin[t] ; 2π c^2 + (2π a b)/Sin[t] (Cos[t]^2 EllipticF[t, k^2] + Sin[t]^2EllipticE[t, k^2])

54.6901

Prolate case: b→c ; surface of revolution (radius c)

2π c^2 + 2π a c ArcSin[(a^2 - c^2)^(1/2)/a]/((a^2 - c^2)^(1/2)/a)

Oblate case: b→a ; surface of revolution (radius a)

2π a^2 + 2π a c ArcSinh[(a^2 - c^2)^(1/2)/c]/((a^2 - c^2)^(1/2)/c)

c→a ; surface of sphere

4π a^2

As a 'pocket calculator' approximation to the surface area of an ellipsoid, that is better than 1.2% (if a > b ≥ c), we may use:

r = ArcCos[c/a]/(1 - c^2/a^2)^(1/2) ;

2 π (c^2 + a b r (1 - (b^2 - c^2)/(6b^2) r^2 (1 - (3 b^2 + 10 c^2)/(56 b^2) r^2)))

54.743

This formula is shorter and has a smaller more symmetric error than the one given in a previous version of this note (2003).
Note that the first two terms correspond to the prolate surface of revolution.

References

1. Wolfram: The Mathematica Book, Wolfram Media, Inc., Fourth Edition, 1999
2. Gradshteyn/Ryzhik: Table of Integrals, Series and Products, Academic Press, Second Printing, 1981

Created by Mathematica  (January 8, 2009)      Go to Download Page