Come utilizzare l'integrazione numerica per calcolare l'area superficiale di un superellissoide?

Aug 19 2020

Sto lavorando in un'applicazione in cui ho bisogno di calcolare l'area superficiale di un superellissoide. Ho letto che non esiste una soluzione in forma chiusa (vedi qui ), quindi sto cercando di calcolarla usando l'integrazione numerica. Il problema è che trovo risultati diversi utilizzando metodi di integrazione numerica diversi, il che mi fa sospettare una sorta di problema di stabilità numerica.

Quindi, le mie domande sono:

  1. C'è un problema numerico con il mio approccio attuale?

  2. Qualcuno sa di un metodo migliore o di miglioramenti che posso apportare per calcolare la superficie di un superellissoide?

Apprezzerei qualsiasi consiglio su come procedere!

Nota: in precedenza avevo pubblicato questa domanda su Mathematics Stack Exchange e non ho ricevuto alcuna risposta. Quindi, sulla base di un suggerimento di mod, ho spostato la domanda su questo sito.


Di seguito i dettagli della mia attuale procedura:

In coordinate cartesiane, abbiamo un'equazione per descrivere un superellissoide in 3D: $$ \left| \frac{x}{r_1} \right|^k + \left| \frac{y}{r_2} \right|^k + \left| \frac{z}{r_3} \right|^k =1 $$

dove $r_1$, $r_2$, e $r_3$ sono le lunghezze dei raggi lungo il $x$, $y$, e $z$assi, rispettivamente. Il parametro$k$definisce la "forma". Se$k=2$, quindi il superellissoide si trasforma in un ellissoide. Come$k \to \infty$, quindi il superellissoide si trasforma in un cuboide. Per un superellissoide generale, gli esponenti di ogni termine possono essere diversi, ma a me interessa solo il caso in cui siano identici.

Dato $r_1$, $r_2$, e $r_3$, Vorrei calcolare l'area della superficie per l'intermedio $k$, dove non si può fare affidamento su formule esistenti per ellissoidi e cuboidi.

Il mio approccio è stato invece quello di utilizzare una rappresentazione parametrica (come mostrato qui su 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} $$

dove le funzioni $c(\alpha, \beta)$ e $s(\alpha, \beta)$ sono definiti come

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

e abbiamo quello $-\pi \leq u < \pi$ e $-\frac{\pi}{2} \leq v < \frac{\pi}{2}$.

A causa della simmetria del superellissoide, possiamo considerare solo la regione in cui $x, y, z \geq 0$, o $0 \leq u, v \leq \frac{\pi}{2}$. Questa regione corrisponde a un ottavo pezzo del superellissoide perché è centrata sull'origine. Quindi, dovremmo semplicemente moltiplicare il nostro risultato in questa regione per$8$ per recuperare la superficie finale.

Con questa semplificazione che $0 \leq u, v \leq \frac{\pi}{2}$, riscriviamo la forma parametrica come:

$$ \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 questa forma parametrica, l'integrale per l'area della superficie è (secondo questo )

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

dove il vettore di posizione $\vec{x} (u, v) = x(u, v) \hat{i} + y(u, v) \hat{j} + z(u, v) \hat{k}$ e il fattore di $8$proviene dall'argomento della simmetria. Valutare l'espressione in Mathematica e semplificare:

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

Sto usando l'espressione di cui sopra per integrare numericamente e trovare l'area della superficie di un superellissoide. Sto testando il semplice caso in cui$r_1=r_2=r_3=1$. In questo caso, abbiamo una sfera unitaria quando$k=2$ con superficie $4 \pi$. Come$k$ diventa grande, quindi la superficie si avvicina $24$. L'area della superficie calcolata per l'intermedio$k$ dovrebbe rientrare in questi limiti.

Sto codificando in R e ho provato a utilizzare le funzioni di integrazione numerica nei pacchetti pracmae cubature. Tra i metodi di integrazione numerica specifici che ho provato con queste funzioni ci sono: quadratura di Gauss-Kronrod, integrazione multidimensionale adattiva (cubatura) e regola di Simpson.

Le diverse implementazioni di integrazione numerica danno risultati estremamente diversi. La maggior parte di loro fornisce risultati troppo piccoli. Alcuni di loro tornano subito NaNper qualsiasi$k>2$. Solo due o tre dei metodi che ho provato (cubatura e qualche variante della quadratura gaussiana) hanno dato risultati ragionevoli ma sono stati un po 'più lenti di quanto avessi sperato. E tutti i metodi falliscono quando$k$ è grande (a partire da circa $k=60$).

Con questi problemi in mente, ci sono problemi numerici alla base di queste discrepanze tra i metodi di integrazione? Esistono modi per risolvere questi problemi? O ancora meglio, esiste un metodo alternativo per calcolare la superficie del superellissoide che eviti questi problemi?

La trama seguente mostra le sfide che ho incontrato con diversi metodi di integrazione numerica. L'asse orizzontale mostra diversi valori di$k$ dove $k=2$ è una forma ellissoide e $k \to \infty$è una forma cuboide. L'asse verticale mostra la superficie data$k$ e una serie di lunghezze di raggio $r_1, r_2, r_3$. In questo caso, le lunghezze degli assi sono il doppio delle lunghezze del raggio. Il grafico mostra l'area della superficie calcolata in funzione di$k$per lo stesso insieme di lunghezze del raggio utilizzando diversi metodi di integrazione numerica. Questi metodi sono implementati in R attraverso i pacchetti pracma(for dblquad) e cubature(per tutti gli altri).

Le due linee verdi orizzontali tratteggiate segnano le aree superficiali dei casi dei membri finali di $k$. Cioè, come$k \to 2$, l'area della superficie dovrebbe convergere alla linea verde inferiore. Tutti i metodi riproducono questo comportamento in$k=2$. Come$k$diventa grande, la superficie dovrebbe convergere alla linea verde superiore. Chiaramente, questo comportamento non è soddisfatto per la maggior parte dei metodi. Il dblquadmetodo dà i risultati più sensati, ma fallisce per i più grandi$k$.

EDIT: Anche l'integrazione numerica funziona e fallisce in modo simile ad altri metodi che utilizzano NIntegrate di Mathematica. Ma il messaggio di errore è più informativo: "Integrazione numerica che converge troppo lentamente; sospetta uno dei seguenti: singolarità, valore dell'integrazione è 0, integrando altamente oscillatorio o WorkingPrecision troppo piccolo."

Risposte

2 Bort Aug 20 2020 at 00:44

Disclaimer, guardo solo il problema con $r_1=r_2=r_3=r=1$. Ma mi aspetto che si possa generalizzare questo approccio per diversi$r_i$.

Suggerisco la seguente mappatura:

Proietta le superfici di un cubo interno sulla superficie del tuo superellissoide. Questo divide la superficie in 6 parti. A causa della simmetria, lo limiterò ora alla mappatura del lato superiore del cubo interno.

Come proiezione, scegliamo la linea che collega l'origine e un punto sulla superficie. L'intersezione della superficie dei cubi sono le coordinate locali$u,v$. Inoltre, lo limiterò a pari$k$, per evitare segni.

Quindi questo dà $$\lambda \left(\begin{array}{c}u\\v\\z\end{array}\right)=x$$ Se usiamo questo nella definizione $$\lambda^k u^k +\lambda^k v^k +\lambda^k z^k =1$$ Otteniamo $\lambda=\left(\frac{1}{u^k+v^k+z^k}\right)^\frac{1}{k}$. Adesso$z$e il dominio di integrazione è ancora indefinito. Qui calcoliamo la proiezione di uno degli angoli sulla tua superficie, con$$\gamma \left(\begin{array}{c}1\\1\\1\end{array}\right)=x$$ otteniamo $$\gamma=\left(\frac{1}{3}\right)^\frac{1}{k}$$. Questo ci dà il dominio di integrazione in$u\in[-\gamma,\gamma]$ e $v\in[-\gamma,\gamma]$ così come $z=\gamma$.

Quindi otteniamo la seguente parametrizzazione in $u,v$ per il lato superiore del tuo superellissoide $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ che sono tutte espressioni di $k$ ovviamente.

Mathematica fornisce come 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}$$

che può essere integrato con k = 100 senza problemi.

Per dispari $k$bisogna controllare attentamente i segni delle espressioni. Questo non dovrebbe essere troppo difficile da risolvere.

Per $k=4$, il rosso mostra una parte del superellissoide che è parametrizzata in u, v. Viene mostrata la metà arancione del modulo completo, il cubo interno e la linea di proiezione.