Как использовать численное интегрирование для вычисления площади поверхности суперэллипсоида?
Я работаю в приложении, в котором мне нужно рассчитать площадь поверхности суперэллипсоида. Я читал, что не существует решения в закрытой форме (см. Здесь ), поэтому я пытаюсь вычислить его, используя численное интегрирование. Проблема в том, что я нахожу разные результаты, используя разные методы численного интегрирования, что заставляет меня подозревать какую-то проблему численной стабильности.
Итак, мои вопросы:
Есть ли числовая проблема с моим нынешним подходом?
Кто-нибудь знает лучший метод или улучшения, которые я могу внести для расчета площади поверхности суперэллипсоида?
Буду признателен за любой совет, как действовать дальше!
Примечание. Я ранее публиковал этот вопрос на сайте Mathematics Stack Exchange и не получил никаких ответов. Поэтому, основываясь на предложении мода, я переместил вопрос на этот сайт.
Подробности моей текущей процедуры ниже:
В декартовых координатах у нас есть уравнение для описания суперэллипсоида в 3D: $$ \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$, то суперэллипсоид превращается в эллипсоид. Так как$k \to \infty$, то суперэллипсоид превращается в кубоид. Для общего суперэллипсоида показатели степени на каждом члене могут быть разными, но меня интересует только случай, в котором они идентичны.
Данный $r_1$, $r_2$, и $r_3$, Я хотел бы вычислить площадь поверхности для промежуточных $k$, где нельзя полагаться на существующие формулы для эллипсоидов и кубоидов.
Мой подход заключался в том, чтобы вместо этого использовать параметрическое представление (как показано здесь, в Википедии):
$$ \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}$.
В силу симметрии суперэллипсоида мы можем рассматривать только область, где $x, y, z \geq 0$, или $0 \leq u, v \leq \frac{\pi}{2}$. Эта область соответствует одной восьмой части суперэллипсоида, поскольку ее центр находится в начале координат. Затем мы просто умножили бы наш результат в этой области на$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, сильно колеблющееся подынтегральное выражение или рабочая точность слишком мала».
Ответы
Отказ от ответственности, я просто смотрю на проблему с $r_1=r_2=r_3=r=1$. Но я ожидаю, что этот подход можно обобщить для разных$r_i$.
Предлагаю следующее отображение:
Спроецируйте поверхности внутреннего куба на поверхность суперэллипсоида. Это делит поверхность на 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$ для верхней стороны вашего суперэллипсоида $$x=\left(\begin{array}{c}\lambda(u,v)u\\\lambda(u,v)v\\\lambda(u,v)\gamma\end{array}\right)$$ которые все выражения $k$ конечно.
Mathematica дает подынтегральное выражение: $$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$, красным показана часть суперэллипсоида, параметризованная по u, v. Показана оранжевая половина полной формы, внутренний куб и линия проекции.
