Bir süperelipsoidin yüzey alanını hesaplamak için sayısal entegrasyon nasıl kullanılır?

Question

Bir süperelipsoidin yüzey alanını hesaplamam gereken bir uygulamada çalışıyorum. Kapalı form çözümü olmadığını okudum ( buraya bakın ), bu yüzden bunu sayısal entegrasyon kullanarak hesaplamaya çalışıyorum. Sorun şu ki, farklı sayısal entegrasyon yöntemlerini kullanarak farklı sonuçlar buluyorum, bu da beni bir tür sayısal kararlılık sorunundan şüpheleniyor.

Yani sorularım:

1. Mevcut yaklaşımımda sayısal bir sorun var mı?

2. Bir süperelipsoidin yüzey alanını hesaplamak için yapabileceğim daha iyi bir yöntem veya iyileştirmeler bilen var mı?

Nasıl ilerleyeceğime dair her türlü tavsiyeyi takdir ediyorum!

Not: Bu soruyu daha önce Matematik Yığın Değişiminde yayınlamıştım ve herhangi bir yanıt almadım. Bu yüzden bir mod önerisine dayanarak soruyu bu siteye taşıdım.

---

Mevcut prosedürümün detayları aşağıda:

Kartezyen koordinatlarda, bir süperelipsoidi 3B olarak tanımlamak için bir denklemimiz var: $$ \left| \frac{x}{r_1} \right|^k + \left| \frac{y}{r_2} \right|^k + \left| \frac{z}{r_3} \right|^k =1 $$

nerede $r_1$, $r_2$, ve $r_3$ boyunca yarıçapların uzunlukları $x$, $y$, ve $z$sırasıyla eksenler. Parametre$k$"şekli" tanımlar. Eğer$k=2$sonra süperelipsoid bir elipsoide dönüşür. Gibi$k \to \infty$sonra süperelipsoid bir küboide dönüşür. Genel bir süperelipsoid için, her terimdeki üsler farklı olabilir, ancak ben sadece özdeş oldukları durumla ilgileniyorum.

Verilen $r_1$, $r_2$, ve $r_3$, Orta seviye için yüzey alanını hesaplamak istiyorum $k$, elipsoidler ve küboidler için mevcut formüllere güvenilemeyecek.

Benim yaklaşımım bunun yerine parametrik bir sunum kullanmak oldu ( burada Wikipedia'da gösterildiği gibi ):

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

fonksiyonlar nerede $c(\alpha, \beta)$ ve $s(\alpha, \beta)$ olarak tanımlanır

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

ve bizde var $-\pi \leq u < \pi$ ve $-\frac{\pi}{2} \leq v < \frac{\pi}{2}$.

Süperelipsoidin simetrisi nedeniyle, yalnızca bulunduğu bölgeyi düşünebiliriz. $x, y, z \geq 0$veya $0 \leq u, v \leq \frac{\pi}{2}$. Bu bölge, süperelipsoidin sekizde bir parçasına karşılık gelir, çünkü orijine merkezlenmiştir. Ardından, bu bölgedeki sonucumuzu basitçe$8$ son yüzey alanını almak için.

Bu sadeleştirme ile $0 \leq u, v \leq \frac{\pi}{2}$, parametrik formu şu şekilde yeniden yazıyoruz:

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

Bu parametrik formu ile yüzey alanı için entegre göre (olup bu )

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

pozisyon vektörü nerede $\vec{x} (u, v) = x(u, v) \hat{i} + y(u, v) \hat{j} + z(u, v) \hat{k}$ ve faktörü $8$simetri argümanından geldi. Mathematica'daki ifadenin değerlendirilmesi ve basitleştirilmesi:

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

Bir süperelipsoidin yüzey alanını sayısal olarak bütünleştirmek ve bulmak için yukarıdaki ifadeyi kullanıyorum. Basit durumu test ediyorum$r_1=r_2=r_3=1$. Bu durumda, bir birim küremiz var$k=2$ yüzey alanı ile $4 \pi$. Gibi$k$ genişler, sonra yüzey alanı yaklaşır $24$. Orta seviye için hesaplanan yüzey alanı$k$ bu sınırlar içinde olmalıdır.

R'de kodluyorum ve pracmave cubaturepaketlerinde sayısal entegrasyon işlevlerini kullanmaya çalıştım . Bu fonksiyonlarla denediğim belirli sayısal entegrasyon yöntemleri arasında: Gauss-Kronrod kuadratürü, uyarlanabilir çok boyutlu entegrasyon (kübik) ve Simpson kuralı vardır.

Farklı sayısal entegrasyon uygulamaları çılgınca farklı sonuçlar verir. Çoğu, çok küçük sonuçlar verir. Bazıları NaNherhangi biri için doğrudan geri döner$k>2$. Denediğim yöntemlerden sadece iki veya üçü (kübik ve bazı Gauss kuadratürü varyantları) mantıklı sonuçlar verdi ama umduğumdan biraz daha yavaş çalışıyor. Ve tüm yöntemler başarısız olduğunda$k$ büyük (yaklaşık olarak $k=60$).

Bu sorunlar göz önünde bulundurularak, entegrasyon yöntemleri arasındaki bu farklılıkların arkasında sayısal sorunlar var mı? Bu sorunları çözebileceğim yollar var mı? Ya da daha iyisi, süperelipsoid yüzey alanını hesaplamak için bu sorunları önleyen alternatif bir yöntem var mı?

Aşağıdaki grafik, farklı sayısal entegrasyon yöntemlerinde karşılaştığım zorlukları göstermektedir. Yatay eksen farklı değerleri gösterir$k$ nerede $k=2$ bir elipsoid şeklidir ve $k \to \infty$küboid bir şekildir. Dikey eksen verilen yüzey alanını gösterir$k$ ve bir dizi yarıçap uzunluğu $r_1, r_2, r_3$. Bu durumda, eksen uzunlukları yarıçap uzunluklarının iki katıdır. Çizim, hesaplanan yüzey alanını bir fonksiyonu olarak gösterir.$k$aynı yarıçap uzunluğu kümesi için çeşitli sayısal entegrasyon yöntemleri kullanılarak. Bu yöntemler R'de pracma(for dblquad) ve cubature(diğerleri için ) paketler aracılığıyla uygulanır .

Kesikli iki yatay yeşil çizgi, uç eleman kasalarının yüzey alanlarını gösterir. $k$. Yani$k \to 2$yüzey alanı alttaki yeşil çizgiye yakınsamalıdır. Tüm yöntemler bu davranışı şurada yeniden üretir:$k=2$. Gibi$k$genişlediğinde, yüzey alanı üstteki yeşil çizgiye yaklaşmalıdır. Açıkça, bu davranış çoğu yöntem için yerine getirilmemiştir. dblquadYöntem en mantıklı sonuç verir, ancak daha büyük başarısız$k$.

DÜZENLEME: Sayısal entegrasyon, Mathematica'nın NIntegrate kullanan diğer yöntemlere benzer şekilde performans gösterir ve başarısız olur. Ancak hata mesajı daha bilgilendirici: "Sayısal entegrasyon çok yavaş yakınsıyor; şunlardan birinden şüphelenin: tekillik, entegrasyonun değeri 0, yüksek salınımlı integrand veya WorkingPrecision çok küçük."

Answer 1

Feragatname, sadece soruna bakıyorum $r_1=r_2=r_3=r=1$. Ancak, bu yaklaşımı farklı$r_i$.

Aşağıdaki eşleştirmeyi öneririm:

Bir iç küpün yüzeylerini süperelipsoidinizin yüzeyine yansıtın. Bu, yüzeyi 6 parçaya böler. Simetri nedeniyle, bunu şimdi iç küpün üst tarafının haritalamasıyla sınırlayacağım.

İzdüşüm olarak, başlangıç ​​noktası ile yüzeydeki bir noktayı birleştiren çizgiyi seçiyoruz. Küp yüzeyinin kesişimi yerel koordinatlardır$u,v$. Dahası, bunu bile sınırlayacağım$k$, işaretlerden kaçınmak için.

Yani bu verir $$\lambda \left(\begin{array}{c}u\\v\\z\end{array}\right)=x$$ Bunu tanımda kullanırsak $$\lambda^k u^k +\lambda^k v^k +\lambda^k z^k =1$$ Elde ederiz $\lambda=\left(\frac{1}{u^k+v^k+z^k}\right)^\frac{1}{k}$. Şimdi$z$ve entegrasyon alanı hala tanımlanmamıştır. Burada, köşelerden birinin yüzeyinize izdüşümünü hesaplıyoruz.$$\gamma \left(\begin{array}{c}1\\1\\1\end{array}\right)=x$$ elde ederiz $$\gamma=\left(\frac{1}{3}\right)^\frac{1}{k}$$. Bu bize entegrasyon alanını verir$u\in[-\gamma,\gamma]$ ve $v\in[-\gamma,\gamma]$ Hem de $z=\gamma$.

Bu nedenle aşağıdaki parametrelendirmeyi elde ederiz $u,v$ süperelipsoidinizin üst tarafı için $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ bunların hepsi ifadeleridir $k$ elbette.

Mathematica integrand olarak verir: $$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}$$

k = 100 ile sorunsuz bir şekilde entegre edilebilir.

Garip için $k$İfadelerin işaretleri dikkatlice kontrol edilmelidir. Bunu düzeltmek çok zor olmamalı.

İçin $k=4$kırmızı, süperelipsoidin u, v'de parametrik hale getirilmiş bir bölümünü gösterir. Tam formun turuncu yarısı, iç küp ve izdüşüm çizgisi gösterilir.