¿Cómo utilizar la integración numérica para calcular el área de superficie de un superelipsoide?
Estoy trabajando en una aplicación en la que necesito calcular el área de superficie de un superelipsoide. He leído que no hay una solución de forma cerrada (ver aquí ), así que estoy tratando de calcularla usando la integración numérica. El problema es que estoy encontrando diferentes resultados utilizando diferentes métodos de integración numérica, lo que me hace sospechar algún tipo de problema de estabilidad numérica.
Entonces, mis preguntas son:
¿Existe un problema numérico con mi enfoque actual?
¿Alguien sabe de un método mejor o de mejoras que pueda realizar para calcular el área de superficie de un superelipsoide?
¡Agradecería cualquier consejo sobre cómo proceder!
Nota: Anteriormente había publicado esta pregunta en Mathematics Stack Exchange y no recibí ninguna respuesta. Entonces, en base a una sugerencia de modificación, moví la pregunta a este sitio.
Detalles de mi procedimiento actual a continuación:
En coordenadas cartesianas, tenemos una ecuación para describir un superelipsoide en 3D: $$ \left| \frac{x}{r_1} \right|^k + \left| \frac{y}{r_2} \right|^k + \left| \frac{z}{r_3} \right|^k =1 $$
dónde $r_1$, $r_2$y $r_3$ son las longitudes de los radios a lo largo del $x$, $y$y $z$ejes, respectivamente. El parámetro$k$define la "forma". Si$k=2$, entonces el superelipsoide se convierte en un elipsoide. Como$k \to \infty$, entonces el superelipsoide se convierte en un cuboide. Para un superelipsoide general, los exponentes de cada término pueden ser diferentes, pero solo me interesa el caso en el que sean idénticos.
Dado $r_1$, $r_2$y $r_3$, Me gustaría calcular el área de superficie para intermedio $k$, donde no se puede confiar en fórmulas existentes para elipsoides y cuboides.
Mi enfoque ha sido, en cambio, utilizar una representación paramétrica (como se muestra aquí en 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} $$
donde las funciones $c(\alpha, \beta)$ y $s(\alpha, \beta)$ se definen como
$$ \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} $$
y tenemos eso $-\pi \leq u < \pi$ y $-\frac{\pi}{2} \leq v < \frac{\pi}{2}$.
Debido a la simetría del superelipsoide, podemos considerar considerar solo la región donde $x, y, z \geq 0$o $0 \leq u, v \leq \frac{\pi}{2}$. Esta región corresponde a una octava parte del superelipsoide porque está centrada en el origen. Entonces, simplemente multiplicaríamos nuestro resultado en esta región por$8$ para recuperar la superficie final.
Con esta simplificación que $0 \leq u, v \leq \frac{\pi}{2}$, reescribimos la forma paramétrica como:
$$ \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} $$
Con esta forma paramétrica, la integral para el área de la superficie es (según esto )
$$ 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 $$
donde el vector de posición $\vec{x} (u, v) = x(u, v) \hat{i} + y(u, v) \hat{j} + z(u, v) \hat{k}$ y el factor de $8$proviene del argumento de la simetría. Evaluar la expresión en Mathematica y simplificar:
$$ 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 $$
Estoy usando la expresión anterior para integrar numéricamente y encontrar el área de superficie de un superelipsoide. Estoy probando el caso simple en el que$r_1=r_2=r_3=1$. En este caso, tenemos una esfera unitaria cuando$k=2$ con superficie $4 \pi$. Como$k$ se vuelve grande, entonces la superficie se acerca $24$. El área de superficie calculada para intermedio$k$ debe estar dentro de esos límites.
Estoy codificando en R y he intentado usar funciones de integración numérica en los paquetes pracma
y cubature
. Entre los métodos de integración numérica específicos que he probado con estas funciones están: cuadratura de Gauss-Kronrod, integración multidimensional adaptativa (cubatura) y la regla de Simpson.
Las diferentes implementaciones de integración numérica dan resultados tremendamente diferentes. La mayoría de ellos dan resultados demasiado pequeños. Algunos de ellos regresan directamente NaN
por cualquier$k>2$. Solo dos o tres de los métodos que probé (cubicación y alguna variante de cuadratura gaussiana) dieron resultados razonables, pero se ejecutaron un poco más lento de lo que esperaba. Y todos los métodos fallan cuando$k$ es grande (a partir de aproximadamente $k=60$).
Con estas cuestiones en mente, ¿hay cuestiones numéricas detrás de estas discrepancias entre los métodos de integración? ¿Existen formas de resolver estos problemas? O mejor aún, ¿existe un método alternativo para calcular el área de superficie del superelipsoide que evite estos problemas?
La siguiente gráfica muestra los desafíos que he encontrado con diferentes métodos de integración numérica. El eje horizontal muestra diferentes valores de$k$ dónde $k=2$ es una forma elipsoide y $k \to \infty$es una forma cuboide. El eje vertical muestra el área de superficie dada$k$ y un conjunto de longitudes de radio $r_1, r_2, r_3$. En este caso, las longitudes de los ejes son el doble de las longitudes de los radios. La gráfica muestra el área de la superficie calculada en función de$k$para el mismo conjunto de longitudes de radio utilizando varios métodos de integración numérica. Estos métodos se implementan en R a través de los paquetes pracma
(para dblquad
) y cubature
(para todos los demás).
Las dos lneas verdes horizontales discontinuas marcan las reas de superficie de los $k$. Es decir, como$k \to 2$, el área de la superficie debe converger con la línea verde inferior. Todos los métodos reproducen este comportamiento en$k=2$. Como$k$se vuelve grande, el área de la superficie debe converger hacia la línea verde superior. Claramente, este comportamiento no se cumple con la mayoría de los métodos. El dblquad
método da los resultados más sensibles, pero falla para grandes$k$.

EDITAR: La integración numérica también funciona y falla de manera similar a otros métodos que utilizan NIntegrate de Mathematica. Pero el mensaje de error es más informativo: "Integración numérica convergiendo demasiado lentamente; sospeche uno de los siguientes: singularidad, el valor de la integración es 0, integrando muy oscilatorio o Precisión de trabajo demasiado pequeña".
Respuestas
Descargo de responsabilidad, solo veo el problema con $r_1=r_2=r_3=r=1$. Pero espero que se pueda generalizar este enfoque para diferentes$r_i$.
Sugiero el siguiente mapeo:
Proyecte las superficies de un cubo interior sobre la superficie de su superelipsoide. Esto divide la superficie en 6 partes. Debido a la simetría, restringiré esto ahora al mapeo del lado superior del cubo interior.
Como proyección, elegimos la línea que conecta el origen y un punto en la superficie. La intersección de la superficie de los cubos son las coordenadas locales$u,v$. Además, restringiré esto incluso$k$, para evitar señales.
Entonces esto da $$\lambda \left(\begin{array}{c}u\\v\\z\end{array}\right)=x$$ Si usamos esto en la definición $$\lambda^k u^k +\lambda^k v^k +\lambda^k z^k =1$$ Obtenemos $\lambda=\left(\frac{1}{u^k+v^k+z^k}\right)^\frac{1}{k}$. Ahora$z$y el dominio de integración aún no está definido. Aquí calculamos la proyección de una de las esquinas sobre su superficie, con$$\gamma \left(\begin{array}{c}1\\1\\1\end{array}\right)=x$$ obtenemos $$\gamma=\left(\frac{1}{3}\right)^\frac{1}{k}$$. Esto nos da el dominio de integración en$u\in[-\gamma,\gamma]$ y $v\in[-\gamma,\gamma]$ tanto como $z=\gamma$.
Por tanto obtenemos la siguiente parametrización en $u,v$ para la parte superior de su superelipsoide $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ que son todas expresiones de $k$ por supuesto.
Mathematica da como integrando: $$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}$$
que se puede integrar con k = 100 sin ningún problema.
Por extraño $k$hay que comprobar cuidadosamente los signos de las expresiones. Esto no debería ser demasiado difícil de solucionar.

Xa $k=4$, el rojo muestra una parte del superelipsoide que está parametrizado en u, v. Se muestra la mitad naranja de la forma completa, el cubo interior y la línea de proyección.
