Comment utiliser l'intégration numérique pour calculer la surface d'un superellipsoïde?

Question

Je travaille dans une application dans laquelle j'ai besoin de calculer la surface d'un superellipsoïde. J'ai lu qu'il n'y avait pas de solution de forme fermée (voir ici ), donc j'essaye de la calculer en utilisant l'intégration numérique. Le problème est que je trouve des résultats différents en utilisant différentes méthodes d'intégration numérique, ce qui me fait suspecter une sorte de problème de stabilité numérique.

Donc, mes questions sont:

1. Y a-t-il un problème numérique avec mon approche actuelle?

2. Est-ce que quelqu'un connaît une meilleure méthode ou des améliorations que je peux apporter pour calculer la surface d'un superellipsoïde?

J'apprécierais tout conseil sur la façon de procéder!

Remarque: j'avais déjà publié cette question sur Mathematics Stack Exchange et je n'ai reçu aucune réponse. Donc, sur la base d'une suggestion de mod, j'ai déplacé la question sur ce site.

---

Détails de ma procédure actuelle ci-dessous:

En coordonnées cartésiennes, nous avons une équation pour décrire un superellipsoïde en 3D: $$ \left| \frac{x}{r_1} \right|^k + \left| \frac{y}{r_2} \right|^k + \left| \frac{z}{r_3} \right|^k =1 $$

où $r_1$, $r_2$, et $r_3$ sont les longueurs des rayons le long du $x$, $y$, et $z$axes, respectivement. Le paramètre$k$définit la «forme». Si$k=2$, puis le superellipsoïde se transforme en ellipsoïde. Comme$k \to \infty$, puis le superellipsoïde se transforme en un cuboïde. Pour un superellipsoïde général, les exposants sur chaque terme peuvent être différents, mais je ne m'intéresse qu'au cas où ils sont identiques.

Donné $r_1$, $r_2$, et $r_3$, Je voudrais calculer la surface pour intermédiaire $k$, où l'on ne peut pas se fier aux formules existantes pour les ellipsoïdes et les cuboïdes.

Mon approche a été d'utiliser à la place une représentation paramétrique (comme indiqué ici sur 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} $$

où les fonctions $c(\alpha, \beta)$ et $s(\alpha, \beta)$ sont définis comme

$$ \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} $$

et nous avons ça $-\pi \leq u < \pi$ et $-\frac{\pi}{2} \leq v < \frac{\pi}{2}$.

En raison de la symétrie du superellipsoïde, on ne peut considérer que la région où $x, y, z \geq 0$, ou $0 \leq u, v \leq \frac{\pi}{2}$. Cette région correspond à un huitième morceau du superellipsoïde car elle est centrée sur l'origine. Ensuite, nous multiplierions simplement notre résultat dans cette région par$8$ pour récupérer la surface finale.

Avec cette simplification que $0 \leq u, v \leq \frac{\pi}{2}$, nous réécrivons la forme paramétrique comme:

$$ \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} $$

Avec cette forme paramétrique, l'intégrale de la surface est (selon cela )

$$ 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 $$

où le vecteur de position $\vec{x} (u, v) = x(u, v) \hat{i} + y(u, v) \hat{j} + z(u, v) \hat{k}$ et le facteur de $8$est venu de l'argument de la symétrie. Évaluer l'expression dans Mathematica et simplifier:

$$ 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 $$

J'utilise l'expression ci-dessus pour intégrer numériquement et trouver la surface d'un superellipsoïde. Je teste le cas simple dans lequel$r_1=r_2=r_3=1$. Dans ce cas, nous avons une sphère unitaire lorsque$k=2$ avec superficie $4 \pi$. Comme$k$ devient grande, alors la surface approche $24$. La surface calculée pour intermédiaire$k$ devrait être dans ces limites.

Je suis en train de coder en R et j'ai essayé d'utiliser des fonctions d'intégration numérique dans les packages pracmaet cubature. Parmi les méthodes d'intégration numérique spécifiques que j'ai essayées avec ces fonctions figurent: la quadrature de Gauss-Kronrod, l'intégration multidimensionnelle adaptative (cubature) et la règle de Simpson.

Les différentes implémentations d'intégration numérique donnent des résultats très différents. La plupart d'entre eux donnent des résultats trop faibles. Certains d'entre eux reviennent directement NaNpour tout$k>2$. Seules deux ou trois des méthodes que j'ai essayées (cubature et une variante de la quadrature gaussienne) ont donné des résultats raisonnables, mais fonctionnent un peu plus lentement que je ne l'avais espéré. Et toutes les méthodes échouent quand$k$ est grand (à partir d'environ $k=60$).

Avec ces problèmes à l'esprit, y a-t-il des problèmes numériques derrière ces écarts entre les méthodes d'intégration? Existe-t-il des moyens de résoudre ces problèmes? Ou mieux encore, existe-t-il une méthode alternative pour calculer la surface superellipsoïde qui évite ces problèmes?

Le graphique ci-dessous montre les défis que j'ai rencontrés avec différentes méthodes d'intégration numérique. L'axe horizontal montre différentes valeurs de$k$ où $k=2$ est une forme ellipsoïde et $k \to \infty$est une forme cuboïde. L'axe vertical montre la surface donnée$k$ et un ensemble de longueurs de rayon $r_1, r_2, r_3$. Dans ce cas, les longueurs d'axe sont deux fois les longueurs de rayon. Le graphique montre la surface calculée en fonction de$k$pour le même ensemble de longueurs de rayon en utilisant plusieurs méthodes d'intégration numériques. Ces méthodes sont implémentées dans R via les packages pracma(pour dblquad) et cubature(pour tous les autres).

Les deux lignes vertes horizontales en pointillés marquent les surfaces des caissons d'extrémité $k$. Autrement dit, comme$k \to 2$, la surface doit converger vers la ligne verte inférieure. Toutes les méthodes reproduisent ce comportement à$k=2$. Comme$k$devient grande, la surface doit converger vers la ligne verte supérieure. De toute évidence, ce comportement n'est pas respecté pour la plupart des méthodes. La dblquadméthode donne les résultats les plus raisonnables, mais échoue pour les$k$.

EDIT: L'intégration numérique fonctionne et échoue également de la même manière que les autres méthodes utilisant NIntegrate de Mathematica. Mais le message d'erreur est plus informatif: "L'intégration numérique converge trop lentement; suspectez l'un des éléments suivants: singularité, la valeur de l'intégration est 0, intégrande hautement oscillatoire ou WorkingPrecision trop petit."

Answer 1

Avertissement, je regarde juste le problème avec $r_1=r_2=r_3=r=1$. Mais j'espère que l'on peut généraliser cette approche pour différents$r_i$.

Je suggère la cartographie suivante:

Projetez les surfaces d'un cube intérieur sur la surface de votre superellipsoïde. Cela divise la surface en 6 parties. En raison de la symétrie, je vais maintenant limiter cela au mappage de la face supérieure du cube intérieur.

Comme projection, nous choisissons la ligne reliant l'origine et un point sur la surface. L'intersection de la surface des cubes sont les coordonnées locales$u,v$. De plus, je limiterai cela à même$k$, pour éviter les signes.

Donc cela donne $$\lambda \left(\begin{array}{c}u\\v\\z\end{array}\right)=x$$ Si nous utilisons cela dans la définition $$\lambda^k u^k +\lambda^k v^k +\lambda^k z^k =1$$ On obtient $\lambda=\left(\frac{1}{u^k+v^k+z^k}\right)^\frac{1}{k}$. Maintenant$z$et le domaine d'intégration n'est toujours pas défini. Ici, nous calculons la projection de l'un des coins sur votre surface, avec$$\gamma \left(\begin{array}{c}1\\1\\1\end{array}\right)=x$$ on obtient $$\gamma=\left(\frac{1}{3}\right)^\frac{1}{k}$$. Cela nous donne le domaine d'intégration dans$u\in[-\gamma,\gamma]$ et $v\in[-\gamma,\gamma]$ aussi bien que $z=\gamma$.

On obtient donc la paramétrisation suivante dans $u,v$ pour la face supérieure de votre superellipsoïde $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ qui sont toutes des expressions de $k$ bien sûr.

Mathematica donne comme intégrande: $$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}$$

qui peut être intégré avec k = 100 sans aucun problème.

Pour bizarre $k$il faut vérifier soigneusement les signes des expressions. Cela ne devrait pas être trop difficile à résoudre.

Pour $k=4$, le rouge montre une partie du superellipsoïde qui est paramétrée dans u, v. La moitié orange de la forme complète, le cube intérieur et la ligne de projection sont affichés.