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

Question

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."

Answer 1

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.