จะใช้การรวมเชิงตัวเลขเพื่อคำนวณพื้นที่ผิวของ superellipsoid ได้อย่างไร?
ฉันกำลังทำงานในแอปพลิเคชันที่ฉันต้องการคำนวณพื้นที่ผิวของซูเปอร์เรลลิปสอยด์ ฉันได้อ่านพบว่าไม่มีโซลูชันรูปแบบปิด (ดูที่นี่ ) ดังนั้นฉันจึงพยายามคำนวณโดยใช้การรวมตัวเลข ปัญหาคือฉันพบผลลัพธ์ที่แตกต่างกันโดยใช้วิธีการรวมตัวเลขที่แตกต่างกันซึ่งทำให้ฉันสงสัยว่ามีปัญหาความเสถียรของตัวเลขบางอย่าง
ดังนั้นคำถามของฉันคือ:
มีปัญหาด้านตัวเลขกับแนวทางปัจจุบันของฉันหรือไม่?
มีใครทราบวิธีการหรือการปรับปรุงที่ดีกว่าที่ฉันสามารถทำได้เพื่อคำนวณพื้นที่ผิวของ superellipsoid หรือไม่?
ฉันจะขอบคุณคำแนะนำใด ๆ ในการดำเนินการต่อ!
หมายเหตุ: ก่อนหน้านี้ฉันเคยโพสต์คำถามนี้ใน Mathematics Stack Exchange และไม่ได้รับคำตอบใด ๆ ตามคำแนะนำของ mod ฉันจึงย้ายคำถามไปที่ไซต์นี้
รายละเอียดขั้นตอนปัจจุบันของฉันด้านล่าง:
ในพิกัดคาร์ทีเซียนเรามีสมการเพื่ออธิบายซูเปอร์เรลลิปสอยด์ในรูปแบบ 3 มิติ: $$ \left| \frac{x}{r_1} \right|^k + \left| \frac{y}{r_2} \right|^k + \left| \frac{z}{r_3} \right|^k =1 $$
ที่ไหน $r_1$, $r_2$และ $r_3$ คือความยาวของรัศมีตามแนว $x$, $y$และ $z$แกนตามลำดับ พารามิเตอร์$k$กำหนด "รูปร่าง" ถ้า$k=2$จากนั้น superellipsoid จะกลายเป็นทรงรี เช่น$k \to \infty$จากนั้น superellipsoid จะกลายเป็นทรงลูกบาศก์ สำหรับ superellipsoid ทั่วไปเลขชี้กำลังของแต่ละเทอมอาจแตกต่างกัน แต่ฉันสนใจเฉพาะกรณีที่เหมือนกันเท่านั้น
ให้ $r_1$, $r_2$และ $r_3$ฉันต้องการคำนวณพื้นที่ผิวสำหรับตัวกลาง $k$โดยที่เราไม่สามารถพึ่งพาสูตรที่มีอยู่สำหรับทรงรีและทรงลูกบาศก์ได้
แนวทางของฉันคือใช้การแสดงพารามิเตอร์แทน (ดังที่แสดงไว้ที่นี่ใน 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} $$
ที่ฟังก์ชั่น $c(\alpha, \beta)$ และ $s(\alpha, \beta)$ ถูกกำหนดให้เป็น
$$ \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} $$
และเรามีสิ่งนั้น $-\pi \leq u < \pi$ และ $-\frac{\pi}{2} \leq v < \frac{\pi}{2}$.
เนื่องจากความสมมาตรของ superellipsoid เราสามารถพิจารณาเฉพาะพื้นที่ที่ $x, y, z \geq 0$, หรือ $0 \leq u, v \leq \frac{\pi}{2}$. ภูมิภาคนี้สอดคล้องกับส่วนหนึ่งในแปดของ superellipsoid เนื่องจากมีศูนย์กลางอยู่ที่ต้นกำเนิด จากนั้นเราจะคูณผลลัพธ์ของเราในภูมิภาคนี้ด้วย$8$ เพื่อดึงพื้นที่ผิวสุดท้าย
ด้วยความเรียบง่ายนี้เองที่ $0 \leq u, v \leq \frac{\pi}{2}$เราเขียนรูปแบบพารามิเตอร์ใหม่เป็น:
$$ \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} $$
ด้วยรูปแบบพาราเมตริกนี้อินทิกรัลสำหรับพื้นที่ผิวคือ (ตามนี้ )
$$ 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 $$
ที่เวกเตอร์ตำแหน่ง $\vec{x} (u, v) = x(u, v) \hat{i} + y(u, v) \hat{j} + z(u, v) \hat{k}$ และปัจจัยของ $8$มาจากอาร์กิวเมนต์สมมาตร การประเมินนิพจน์ใน Mathematica และทำให้เข้าใจง่าย:
$$ 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 $$
ฉันกำลังใช้นิพจน์ข้างต้นเพื่อรวมตัวเลขและค้นหาพื้นที่ผิวของซูเปอร์เรลลิปสอยด์ ฉันกำลังทดสอบกรณีง่ายๆซึ่ง$r_1=r_2=r_3=1$. ในกรณีนี้เรามีหน่วยทรงกลมเมื่อ$k=2$ ด้วยพื้นที่ผิว $4 \pi$. เช่น$k$ กลายเป็นขนาดใหญ่จากนั้นพื้นที่ผิวจะเข้าใกล้ $24$. พื้นที่ผิวที่คำนวณได้สำหรับตัวกลาง$k$ ควรอยู่ในขอบเขตเหล่านั้น
ฉันกำลังเข้ารหัสใน R และได้พยายามใช้ฟังก์ชันการรวมตัวเลขในpracmaและcubatureแพ็คเกจ ในบรรดาวิธีการรวมเชิงตัวเลขที่เฉพาะเจาะจงที่ฉันได้ลองใช้กับฟังก์ชันเหล่านี้ ได้แก่ การสร้างพื้นที่สี่เหลี่ยมจัตุรัสแบบเกาส์ - โครนรอดการรวมหลายมิติแบบปรับได้ (คิวบิก) และกฎของซิมป์สัน
การใช้งานการรวมตัวเลขที่แตกต่างกันให้ผลลัพธ์ที่แตกต่างกันอย่างมาก ส่วนใหญ่ให้ผลที่เล็กเกินไป บางส่วนของพวกเขาขึ้นตรงกลับNaNสำหรับการใด ๆ$k>2$. มีเพียงสองหรือสามวิธีที่ฉันลองใช้ (รูปสี่เหลี่ยมจัตุรัสและรูปสี่เหลี่ยมแบบเกาส์เซียนบางรูปแบบ) ให้ผลลัพธ์ที่สมเหตุสมผล แต่ทำงานช้ากว่าที่ฉันคาดไว้เล็กน้อย และทุกวิธีล้มเหลวเมื่อ$k$ มีขนาดใหญ่ (เริ่มจากประมาณ $k=60$).
เมื่อคำนึงถึงประเด็นเหล่านี้แล้วมีปัญหาด้านตัวเลขที่อยู่เบื้องหลังความคลาดเคลื่อนระหว่างวิธีการผสานรวมหรือไม่ มีวิธีใดบ้างที่ฉันสามารถแก้ไขปัญหาเหล่านี้ได้ หรือดีกว่านั้นมีวิธีอื่นในการคำนวณพื้นที่ผิวเหนือดินที่หลีกเลี่ยงปัญหาเหล่านี้หรือไม่?
พล็อตด้านล่างแสดงให้เห็นถึงความท้าทายที่ฉันได้พบกับวิธีการรวมตัวเลขที่แตกต่างกัน แกนนอนแสดงค่าต่างๆของ$k$ ที่ไหน $k=2$ เป็นรูปทรงรีและ $k \to \infty$เป็นรูปทรงลูกบาศก์ แกนแนวตั้งแสดงพื้นที่ผิวที่กำหนด$k$ และชุดของความยาวรัศมี $r_1, r_2, r_3$. ในกรณีนี้ความยาวแกนจะเป็นสองเท่าของความยาวรัศมี พล็อตแสดงพื้นที่ผิวที่คำนวณตามฟังก์ชันของ$k$สำหรับความยาวรัศมีชุดเดียวกันโดยใช้วิธีการรวมตัวเลขหลายวิธี วิธีการเหล่านี้ถูกนำไปใช้ใน R ผ่านแพ็คเกจpracma(สำหรับdblquad) และcubature(สำหรับอื่น ๆ ทั้งหมด)
เส้นสีเขียวแนวนอนสองเส้นประทำเครื่องหมายพื้นที่ผิวของกรณีสมาชิกปลายทางของ $k$. นั่นคือเป็น$k \to 2$พื้นที่ผิวควรมาบรรจบกับเส้นสีเขียวด้านล่าง วิธีการทั้งหมดทำให้เกิดพฤติกรรมนี้ที่$k=2$. เช่น$k$กลายเป็นขนาดใหญ่พื้นที่ผิวควรมาบรรจบกับเส้นสีเขียวด้านบน เห็นได้ชัดว่าพฤติกรรมนี้ใช้ไม่ได้กับวิธีการส่วนใหญ่ dblquadวิธีให้ผลที่เหมาะสมที่สุด แต่ล้มเหลวขนาดใหญ่$k$.
แก้ไข: การรวมเชิงตัวเลขยังดำเนินการและล้มเหลวเช่นเดียวกันกับวิธีการอื่น ๆ โดยใช้ NIntegrate ของ Mathematica แต่ข้อความแสดงข้อผิดพลาดนั้นให้ข้อมูลมากกว่า: "การรวมเชิงตัวเลขมาบรรจบกันช้าเกินไปสงสัยอย่างใดอย่างหนึ่งต่อไปนี้: ความเป็นเอกฐานค่าของการรวมเป็น 0 ค่าปริพันธ์ที่มีความผันผวนสูงหรือ WorkingPrecision น้อยเกินไป"
คำตอบ
ข้อจำกัดความรับผิดชอบฉันเพิ่งดูปัญหาด้วย $r_1=r_2=r_3=r=1$. แต่ฉันคาดหวังว่าเราสามารถสรุปแนวทางนี้ให้แตกต่างกันได้$r_i$.
ฉันขอแนะนำการทำแผนที่ต่อไปนี้:
ฉายพื้นผิวของลูกบาศก์ภายในลงบนพื้นผิวของ Superellipsoide ของคุณ สิ่งนี้แบ่งพื้นผิวออกเป็น 6 ส่วน เนื่องจากความสมมาตรตอนนี้ฉันจะ จำกัด สิ่งนี้ไว้ที่การแมปด้านบนของคิวบ์ภายใน
ในการฉายภาพเราเลือกเส้นที่เชื่อมต่อจุดเริ่มต้นและจุดบนพื้นผิว จุดตัดของพื้นผิวลูกบาศก์คือพิกัดท้องถิ่น$u,v$. นอกจากนี้ฉันจะ จำกัด สิ่งนี้ไว้ที่คู่$k$เพื่อหลีกเลี่ยงสัญญาณ
ดังนั้นสิ่งนี้จะช่วยให้ $$\lambda \left(\begin{array}{c}u\\v\\z\end{array}\right)=x$$ ถ้าเราใช้คำจำกัดความนี้ $$\lambda^k u^k +\lambda^k v^k +\lambda^k z^k =1$$ เราได้รับ $\lambda=\left(\frac{1}{u^k+v^k+z^k}\right)^\frac{1}{k}$. ตอนนี้$z$และโดเมนการรวมยังไม่ได้กำหนด ที่นี่เราคำนวณการฉายภาพของมุมใดมุมหนึ่งบนพื้นผิวของคุณด้วย$$\gamma \left(\begin{array}{c}1\\1\\1\end{array}\right)=x$$ เราได้รับ $$\gamma=\left(\frac{1}{3}\right)^\frac{1}{k}$$. สิ่งนี้ทำให้เรามีโดเมนการรวมใน$u\in[-\gamma,\gamma]$ และ $v\in[-\gamma,\gamma]$ เช่นเดียวกับ $z=\gamma$.
ดังนั้นเราจึงได้รับพารามิเตอร์ต่อไปนี้ใน $u,v$ สำหรับด้านบนของ superellipsoid ของคุณ $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ ซึ่งเป็นนิพจน์ทั้งหมดของ $k$ แน่นอน.
Mathematica ให้เป็น 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}$$
ซึ่งสามารถรวมเข้ากับ k = 100 ได้โดยไม่มีปัญหาใด ๆ
สำหรับคี่ $k$เราต้องตรวจสอบสัญญาณของการแสดงออกอย่างรอบคอบ สิ่งนี้ไม่ควรแก้ไขยากเกินไป
สำหรับ $k=4$, สีแดงแสดงส่วนหนึ่งของ superellipsoid ซึ่งเป็น parametrized ใน u, v. ครึ่งสีส้มของรูปแบบเต็มลูกบาศก์ภายในและเส้นฉายจะปรากฏขึ้น
