Bagaimana cara menggunakan integrasi numerik untuk menghitung luas permukaan superellipsoid?
Saya bekerja dalam aplikasi di mana saya perlu menghitung luas permukaan superellipsoid. Saya telah membaca bahwa tidak ada solusi bentuk tertutup (lihat di sini ), jadi saya mencoba menghitungnya menggunakan integrasi numerik. Masalahnya adalah saya menemukan hasil yang berbeda menggunakan metode integrasi numerik yang berbeda, yang membuat saya mencurigai beberapa jenis masalah stabilitas numerik.
Jadi, pertanyaan saya adalah:
Apakah ada masalah numerik dengan pendekatan saya saat ini?
Adakah yang tahu tentang metode atau perbaikan yang lebih baik yang dapat saya lakukan untuk menghitung luas permukaan superellipsoid?
Saya sangat menghargai saran apa pun tentang cara melanjutkan!
Catatan: Saya sebelumnya telah memposting pertanyaan ini di Mathematics Stack Exchange, dan tidak menerima tanggapan apa pun. Jadi berdasarkan saran mod, saya memindahkan pertanyaan ke situs ini.
Rincian prosedur saya saat ini di bawah:
Dalam koordinat Kartesius, kami memiliki persamaan untuk mendeskripsikan superellipsoid dalam 3D: $$ \left| \frac{x}{r_1} \right|^k + \left| \frac{y}{r_2} \right|^k + \left| \frac{z}{r_3} \right|^k =1 $$
dimana $r_1$, $r_2$, dan $r_3$ adalah panjang jari-jari di sepanjang $x$, $y$, dan $z$sumbu, masing-masing. Parameternya$k$mendefinisikan "bentuk". Jika$k=2$, kemudian superellipsoid berubah menjadi ellipsoid. Sebagai$k \to \infty$, kemudian superellipsoid berubah menjadi berbentuk kubus. Untuk superellipsoid umum, eksponen pada setiap suku bisa berbeda, tetapi saya hanya tertarik pada kasus di mana mereka identik.
Diberikan $r_1$, $r_2$, dan $r_3$, Saya ingin menghitung luas permukaan untuk menengah $k$, di mana seseorang tidak dapat mengandalkan rumus yang ada untuk ellipsoid dan kuboid.
Pendekatan saya adalah menggunakan representasi parametrik (seperti yang ditunjukkan di sini di 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} $$
dimana fungsinya $c(\alpha, \beta)$ dan $s(\alpha, \beta)$ didefinisikan sebagai
$$ \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} $$
dan kami memilikinya $-\pi \leq u < \pi$ dan $-\frac{\pi}{2} \leq v < \frac{\pi}{2}$.
Karena kesimetrisan superellipsoid, kita hanya dapat mempertimbangkan wilayah di mana $x, y, z \geq 0$, atau $0 \leq u, v \leq \frac{\pi}{2}$. Wilayah ini sesuai dengan seperdelapan bagian superellipsoid karena berpusat pada asalnya. Kemudian, kami hanya akan mengalikan hasil kami di wilayah ini dengan$8$ untuk mendapatkan kembali luas permukaan akhir.
Dengan penyederhanaan itu $0 \leq u, v \leq \frac{\pi}{2}$, kami menulis ulang bentuk parametrik sebagai:
$$ \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} $$
Dengan bentuk parametrik ini, integral untuk luas permukaan adalah (menurut ini )
$$ 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 $$
dimana posisi vektor $\vec{x} (u, v) = x(u, v) \hat{i} + y(u, v) \hat{j} + z(u, v) \hat{k}$ dan faktor $8$berasal dari argumen simetri. Mengevaluasi ekspresi dalam Mathematica dan menyederhanakan:
$$ 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 $$
Saya menggunakan ekspresi di atas untuk mengintegrasikan secara numerik dan menemukan luas permukaan superellipsoid. Saya menguji kasus sederhana di mana$r_1=r_2=r_3=1$. Dalam hal ini, kami memiliki bola satuan when$k=2$ dengan luas permukaan $4 \pi$. Sebagai$k$ menjadi besar, kemudian luas permukaan mendekat $24$. Luas permukaan yang dihitung untuk menengah$k$ harus dalam batas tersebut.
Saya membuat kode di R, dan telah mencoba menggunakan fungsi integrasi numerik dalam paket pracmadan cubature. Di antara metode integrasi numerik tertentu yang telah saya coba dengan fungsi ini adalah: kuadrat Gauss-Kronrod, integrasi multidimensi adaptif (kubatur), dan aturan Simpson.
Implementasi integrasi numerik yang berbeda memberikan hasil yang sangat berbeda. Kebanyakan dari mereka memberikan hasil yang terlalu kecil. Beberapa dari mereka langsung kembali NaNuntuk apa pun$k>2$. Hanya dua atau tiga metode yang saya coba (kubatur dan beberapa varian kuadratur Gaussian) yang memberikan hasil yang masuk akal tetapi berjalan sedikit lebih lambat dari yang saya harapkan. Dan semua metode gagal saat$k$ besar (mulai dari sekitar $k=60$).
Dengan mempertimbangkan masalah ini, apakah ada masalah numerik di balik perbedaan antara metode integrasi ini? Adakah cara agar saya dapat menyelesaikan masalah ini? Atau lebih baik lagi, adakah metode alternatif untuk menghitung luas permukaan superellipsoid yang menghindari masalah ini?
Plot di bawah ini menunjukkan tantangan yang saya temui dengan metode integrasi numerik yang berbeda. Sumbu horizontal menunjukkan nilai yang berbeda dari$k$ dimana $k=2$ adalah bentuk elipsoid dan $k \to \infty$adalah bentuk kubus. Sumbu vertikal menunjukkan luas permukaan yang diberikan$k$ dan satu set panjang radius $r_1, r_2, r_3$. Dalam kasus ini, panjang sumbu adalah dua kali panjang jari-jari. Plot menunjukkan luas permukaan yang dihitung sebagai fungsi dari$k$untuk kumpulan panjang radius yang sama menggunakan beberapa metode integrasi numerik. Metode ini diimplementasikan di R melalui paket pracma(untuk dblquad) dan cubature(untuk yang lainnya).
Dua garis hijau horizontal putus-putus menandai luas permukaan kotak anggota-akhir $k$. Artinya, sebagai$k \to 2$, luas permukaan harus bertemu dengan garis hijau bawah. Semua metode mereproduksi perilaku ini di$k=2$. Sebagai$k$menjadi besar, luas permukaan harus menyatu dengan garis hijau atas. Jelas, perilaku ini tidak dipenuhi untuk sebagian besar metode. The dblquadmetode memberikan hasil yang paling masuk akal, tapi gagal untuk yang lebih besar$k$.
EDIT: Integrasi numerik juga melakukan dan gagal serupa dengan metode lain menggunakan NIntegrate Mathematica. Namun pesan kesalahannya lebih informatif: "Integrasi numerik berkumpul terlalu lambat; curigai salah satu dari yang berikut: singularitas, nilai integrasi adalah 0, integrand sangat berosilasi, atau WorkingPrecision terlalu kecil."
Jawaban
Penafian, saya hanya melihat masalah dengan $r_1=r_2=r_3=r=1$. Tetapi saya berharap, seseorang dapat menggeneralisasi pendekatan ini untuk hal yang berbeda$r_i$.
Saya menyarankan pemetaan berikut:
Proyeksikan permukaan kubus interior ke permukaan superellipsoide Anda. Ini membagi permukaan menjadi 6 bagian. Karena kesimetrisannya, saya akan membatasi ini sekarang pada pemetaan sisi atas kubus interior.
Sebagai proyeksi, kami memilih garis yang menghubungkan titik asal dan titik di permukaan. Perpotongan permukaan kubus adalah koordinat lokalnya$u,v$. Selanjutnya, saya akan membatasi ini bahkan$k$, untuk menghindari tanda-tanda.
Jadi ini memberi $$\lambda \left(\begin{array}{c}u\\v\\z\end{array}\right)=x$$ Jika kita menggunakan ini dalam definisi $$\lambda^k u^k +\lambda^k v^k +\lambda^k z^k =1$$ Kami mendapatkan $\lambda=\left(\frac{1}{u^k+v^k+z^k}\right)^\frac{1}{k}$. Sekarang$z$dan domain integrasi masih belum ditentukan. Di sini kami menghitung proyeksi salah satu sudut ke permukaan Anda, dengan$$\gamma \left(\begin{array}{c}1\\1\\1\end{array}\right)=x$$ kami dapatkan $$\gamma=\left(\frac{1}{3}\right)^\frac{1}{k}$$. Ini memberi kita domain integrasi dalam$u\in[-\gamma,\gamma]$ dan $v\in[-\gamma,\gamma]$ sebaik $z=\gamma$.
Karenanya kami mendapatkan parametrization berikut di $u,v$ untuk sisi atas superellipsoid Anda $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ yang semuanya merupakan ekspresi $k$ tentu saja.
Mathematica memberi sebagai 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}$$
yang dapat diintegrasikan dengan k = 100 tanpa masalah.
Aneh $k$seseorang harus hati-hati memeriksa tanda-tanda ekspresi. Ini seharusnya tidak terlalu sulit untuk diperbaiki.
Untuk $k=4$, merah menunjukkan bagian dari superellipsoid yang diparameterisasi dalam u, v. Setengah oranye dari bentuk penuh, kubus interior dan garis proyeksi ditampilkan.
