# How to use numerical integration to calculate the surface area of a superellipsoid?

I am working in an application in which I need to calculate the surface area of a superellipsoid. I have read that there is no closed form solution (see here), so I am trying to compute it using numerical integration. The problem is that I am finding different results using different numerical integration methods, which makes me suspect some kind of numerical stability issue.

So, my questions are:

1. Is there a numerical issue with my current approach?

2. Does anyone know of a better method or improvements I can make to calculate the surface area of a superellipsoid?

I would appreciate any advice on how to proceed!

Note: I had previously posted this question on the Mathematics Stack Exchange, and didn't receive any responses. So based on a mod suggestion, I moved the question to this site.

Details of my current procedure below:

In Cartesian coordinates, we have an equation to describe a superellipsoid in 3D: $$\left| \frac{x}{r_1} \right|^k + \left| \frac{y}{r_2} \right|^k + \left| \frac{z}{r_3} \right|^k =1$$

where $$r_1$$, $$r_2$$, and $$r_3$$ are the lengths of the radii along the $$x$$, $$y$$, and $$z$$ axes, respectively. The parameter $$k$$ defines the "shape." If $$k=2$$, then the superellipsoid turns into an ellipsoid. As $$k \to \infty$$, then the superellipsoid turns into a cuboid. For a general superellipsoid, the exponents on each term can be different, but I am only interested in the case in which they are identical.

Given $$r_1$$, $$r_2$$, and $$r_3$$, I would like to compute surface area for intermediate $$k$$, where one cannot rely on existing formulae for ellipsoids and cuboids.

My approach has been to instead use a parametric representation (as shown here on Wikipedia):

\begin{align} x(u, v)&=r_1 c \left(v, \frac{2}{k} \right) c \left( u, \frac{2}{k} \right) \\ y(u, v)&=r_2 c \left(v, \frac{2}{k} \right) s \left( u, \frac{2}{k} \right) \\ z(u, v)&=r_3 s \left( v, \frac{2}{k} \right) \end{align}

where the functions $$c(\alpha, \beta)$$ and $$s(\alpha, \beta)$$ are defined as

\begin{align} c(\alpha, \beta)&=\mathrm{sgn}(\cos{\alpha}) \left| \cos{\alpha} \right|^\beta \\ s(\alpha, \beta)&=\mathrm{sgn}(\sin{\alpha}) \left| \sin{\alpha} \right|^\beta \end{align}

and we have that $$-\pi \leq u < \pi$$ and $$-\frac{\pi}{2} \leq v < \frac{\pi}{2}$$.

Due to the symmetry of the superellipsoid, we can consider consider only the region where $$x, y, z \geq 0$$, or $$0 \leq u, v \leq \frac{\pi}{2}$$. This region corresponds to a one-eighth chunk of the superellipsoid because it is centered on the origin. Then, we would simply multiply our result in this region by $$8$$ to retrieve the final surface area.

With this simplification that $$0 \leq u, v \leq \frac{\pi}{2}$$, we rewrite the parametric form as:

\begin{align} x(u, v)&=r_1 (\cos{v} \cos{u})^\frac{2}{k} \\ y(u, v)&=r_2 (\cos{v} \sin{u})^\frac{2}{k} \\ z(u, v)&=r_3 (\sin{v})^\frac{2}{k} \end{align}

With this parametric form, the integral for the surface area is (according to this)

$$A=\int \int_S \mathrm{d}S = 8 \int_0^\frac{\pi}{2} \int_0^\frac{\pi}{2} \left| \left| \frac{\partial \vec{x} (u, v)}{\partial u} \times \frac{\partial \vec{x} (u, v)}{\partial v} \right| \right| \mathrm{d}u \ \mathrm{d}v$$

where the position vector $$\vec{x} (u, v) = x(u, v) \hat{i} + y(u, v) \hat{j} + z(u, v) \hat{k}$$ and the factor of $$8$$ has come from the symmetry argument. Evaluating the expression in Mathematica and simplifying:

$$A=\frac{32}{k^2} \int_0^\frac{\pi}{2} \int_0^\frac{\pi}{2} \sqrt{\left(r_2 r_3 \cos{u} (\sin{u} \sin{v} \cos{v})^{\frac{2}{k}-1} \cos^2{v} \right)^2 + \left(r_1 r_3 \sin{u} (\cos{u} \sin{v}\cos{v})^{\frac{2}{k}-1} \cos^2{v}) \right)^2 + \left(r_1 r_2 \sin{v} (\sin{u} \cos{u} \cos{v})^{\frac{2}{k}-1} (\cos{v})^\frac{2}{k} \right)^2} \mathrm{d}u \ \mathrm{d}v$$

I am using the above expression to numerically integrate and find the surface area of a superellipsoid. I am testing the simple case in which $$r_1=r_2=r_3=1$$. In this case, we have a unit sphere when $$k=2$$ with surface area $$4 \pi$$. As $$k$$ becomes large, then the surface area approaches $$24$$. The computed surface area for intermediate $$k$$ should be within those limits.

I am coding in R, and have tried to use numerical integration functions in the pracma and cubature packages. Among the specific numerical integration methods I have tried with these functions are: Gauss-Kronrod quadrature, adaptive multidimensional integration (cubature), and Simpson's rule.

The different numerical integration implementations give wildly different results. Most of them give results that are too small. Some of them straight up return NaN for any $$k>2$$. Only two or three of the methods I tried (cubature and some variant of Gaussian quadrature) gave sensible results but run a bit slower than I had hoped. And all of the methods fail when $$k$$ is large (starting from about $$k=60$$).

With these issues in mind, are there numerical issues behind these discrepancies between integration methods? Are there ways I can resolve these issues? Or even better, is there an alternative method to calculate superellipsoid surface area that avoids these issues?

The plot below shows the challenges I have come across with different numerical integration methods. The horizontal axis shows different values of $$k$$ where $$k=2$$ is an ellipsoid shape and $$k \to \infty$$ is a cuboid shape. The vertical axis shows surface area given $$k$$ and a set of radius lengths $$r_1, r_2, r_3$$. In this case, axis lengths are twice the radius lengths. The plot shows computed surface area as a function of $$k$$ for the same set of radius lengths using several numerical integration methods. These methods are implemented in R through the packages pracma (for dblquad) and cubature (for all others).

The two dashed horizontal green lines mark the surface areas of the end-member cases of $$k$$. That is, as $$k \to 2$$, surface area should converge to the bottom green line. All of the methods reproduce this behavior at $$k=2$$. As $$k$$ becomes large, surface area should converge to the top green line. Clearly, this behavior is not fulfilled for most methods. The dblquad method gives the most sensible results, but fails for larger $$k$$. EDIT: Numerical integration also performs and fails similarly to other methods using Mathematica's NIntegrate. But the error message is more informative: "Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small."

• Your expression for the area seems to go to 0 as the value of k increases for constant r (the sines and cosines are bounded in [-1,1]) That doesn't seem right. Are you confident that the expression is correct? – Biswajit Banerjee Aug 18 '20 at 21:23
• @BiswajitBanerjee That is a good point. I went back to Mathematica and checked the integrand, but I got the same result. So I'm confident the integrand is theoretically correct but it does seem strange that it appears it goes to 0 as $k$ increases. I also tried to integrate with Mathematica's NIntegrate. It matches with other methods at $k=2, 3$ but also fails at large $k$. But it gives a more informative error: "Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small." – llorente Aug 19 '20 at 0:33
• I managed to get NIntegrate to work up to $k=20$ with DoupleExponential as method. I suspect, one needs another coordinate mapping than the one you propose for larger $k$. – Bort Aug 19 '20 at 8:30
• The expression you wrote has an integrable singularity in it as either $u$ and $v$ go to 0 for $k > 2$, so it's not surprising that a canned integration routine will fail. For this kind of integrand you'll definitely need adaptive quadrature, but a different parametrization without a singularity would be even better. Another idea -- it's worth calculating the derivative of the answer w.r.t. $k$ around a known value like $k = 2$. This gives you something to benchmark your numerical solution against. – Daniel Shapero Aug 20 '20 at 15:01

## 1 Answer

Disclaimer, I just look at the problem with $$r_1=r_2=r_3=r=1$$. But I expect, that one can generalize this approach for different $$r_i$$.

I suggest the following mapping:

Project the surfaces of an interior cube onto the surface of your superellipsoide. This divides the surface into 6 parts. Because of the symmetry, I will restrict this now to the mapping of the top side of the interior cube.

As the projection, we choose the line connecting the origin and a point on the surface. The intersection of the cubes surface are the local coordinates $$u,v$$. Further, I will restrict this to even $$k$$, to avoid signs.

So this gives $$\lambda \left(\begin{array}{c}u\\v\\z\end{array}\right)=x$$ If we use this in the definition $$\lambda^k u^k +\lambda^k v^k +\lambda^k z^k =1$$ We obtain $$\lambda=\left(\frac{1}{u^k+v^k+z^k}\right)^\frac{1}{k}$$. Now $$z$$ and the integration domain is still undefined. Here we calculate the projection of one of the corners onto your surface, with $$\gamma \left(\begin{array}{c}1\\1\\1\end{array}\right)=x$$ we obtain $$\gamma=\left(\frac{1}{3}\right)^\frac{1}{k}$$. This gives us the integration domain in $$u\in[-\gamma,\gamma]$$ and $$v\in[-\gamma,\gamma]$$ as well as $$z=\gamma$$.

Hence we obtain the following parametrization in $$u,v$$ for the top side of your superellipsoid $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ which are all expressions of $$k$$ of course.

Mathematica gives as integrand: $$3^{-1/k} \sqrt{9^{\frac{1}{k}-1} \left| u^k+v^k+\frac{1}{3}\right| ^{-\frac{2 (k+2)}{k}}+\left| v^{k-1} \left(\frac{1}{u^k+v^k+\frac{1}{3}}\right)^{\frac{k+2}{k}}\right| ^2+\left| u^{k-1} \left(\frac{1}{u^k+v^k+\frac{1}{3}}\right)^{\frac{k+2}{k}}\right| ^2}$$

which can be integrated with k=100 without any problems.

For odd $$k$$ one has to carefully check signs of the expressions. This shouldn't be too hard to fix. For $$k=4$$, red shows a part of the superellipsoid which is parametrized in u,v. Orange half of the full form, the interior cube and the projection line is shown. • That's brilliant. I got the same result for $r_1=r_2=r_3=1$. But I have trouble generalizing to arbitrary $r_i$. I restricted the domain to 1/4 of the cube's top face so that $u, v\geq0$ and we don't have to consider parity of $k$. Then you multiply by 24 to find the full area. I have that $\lambda= \left[ \left(\frac{u}{r_1} \right)^k+\left(\frac{v}{r_2}\right)^k+\left(\frac{z}{r_3}\right)^k)\right]^{-1/k}$ and $\gamma$ is unchanged. Integration region is the same but the lower bounds are 0. When I do this, the area does not match with the expected areas at $k=2$ and $k \to \infty$. – llorente Aug 19 '20 at 23:20
• Well, if you have different $r_i$, you have essentially a cuboid inside instead of a cube. So the resulting full integral is split in the sum of three different integrals. The overall approach for one individual side should still work. – Bort Aug 20 '20 at 8:42
• Thanks for the tip. I have solved the problem now for general $r_i$ and $k$. Yes, like you said, I had to calculate three integrals from the unique combinations of two radius lengths. This corresponds to a "corner" of the superellipsoid, so then you need to multiply by 8 to recover the full surface area. And the integration bounds do actually change, so part of my previous comment was incorrect. – llorente Aug 20 '20 at 21:52