Кратный интеграл по сложной области, как посчитать численно?

[TheEnt]

Чтобы не плодить новый тем, задам вопрос тут.
Допустим есть кратный интеграл $\int_{x_1}^{x_2}dx \int_{y_1(x)}^{y_2(x)}dy f(x,y)$. У меня есть оттестированный алгоритм взятия кратных интегралов с пределами, не зависящими от переменных (прямоугольная область), могу ли я его использовать для взятия интеграла по более сложными функциями, если я сделаю так
$\int_{x_1}^{x_2}dx \int_{y_1(x)}^{y_2(x)}dy f(x,y)=\int_{x_1}^{x_2}dx \int_{y_1}^{y_2}dy \theta(y-y_1(x)) \theta(y_2(x)-y)f(x,y)$,
где $\theta(x)$ - функция Хевисайда, и если известно, что пределы - это достаточно "хорошие" непрерывные функции? Я понимаю, что такой способ взятия интеграла может быть нерациональным, но все же мне интересно, можно ли так делать. На деле у меня есть трех-кратный интеграл, и просто мне не очень хочется писать какой-то сложный метод ради взятия всего одного интеграла...

[Sonic86]

Ну да, все правильно: можно, но нерационально (хотя я не знаю, как рациональнее. Ну есть метод Монте-Карло (хотя вроде даже то же самое) и метод разбиения на области) + еще проблема с вычислением пределов $y_1,y_2$ - нужно искать экстремумы функций $y_j(x)$.

[ewert]

TheEnt в сообщении #509155 писал(а):

Я понимаю, что такой способ взятия интеграла может быть нерациональным, но все же мне интересно, можно ли так делать.

Формально можно, только точность будет крайне низкой -- из-за разрывности подынтегральной функции.

[svv]

TheEnt писал(а):

Чтобы не плодить новый тем, задам вопрос тут.

После этой фразы раза три посмотрел, не является ли эта страничка второй или третьей страницей какой-то старой темы.

АКМ:

svv,

ну не надо по три раза смотреть на мелочи, давайте сразу догадываться.
Ну да, туда надо было добавить, что, мол отделил сюда, сюда добавить, что отделил оттуда, где-то ещё нравоучение автору написать...
Ayez pitié! Вот назначат Вас модератором, посмотрю я на Вас...

[TheEnt]

Спасибо за ответы, я попробовал, точность просто ужасная... Придется делать нормально. Вот только я прочитал такую вещь в книге Калиткина Н.Н. "Численные методы" (стр 112-113), что если кратный интеграл по криволинейной области брать последовательно (т.е. если есть $\int dx \int dy f(x,y)$, сначала берется $F(y)=\int dx f(x,y)$, а потом $\int dy F(y)$), то уже при втором интегрировании интегрирования надо делать разложение с весами полиномов Чебышева. Как поступать с 3 интегралом даже не знаю... Нужно опять раскладывать в ряд с весами по полиномам Чебышева?

-- Пн ноя 28, 2011 17:24:45 --

svv в сообщении #509217 писал(а):

TheEnt писал(а):

Чтобы не плодить новый тем, задам вопрос тут.

После этой фразы раза три посмотрел, не является ли эта страничка второй или третьей страницей какой-то старой темы.

Это была 2 или 3 страничка старой темы, но модераторы выделили мой вопрос в отдельную тему

[TheEnt]

Скажите пожалуйста, что почитать по теме, никак не могу решить задачу. Ещё раз, как звучит вопрос - как численно взять кратный интеграл со сложными пределами? Можно ли использовать метод Монте-Карло? Почти все, что я читал, применимо только для интегрирования по n-мерному параллелепипеду.
Что я читал по моему вопросу - Демидович Б.П., Марон И.А. Основы вычислительной математики - тут по сути описан метод, который я описал в первом посте.
Калиткин Н.Н. - Численные методы - тут описано, как надо брать 2-х кратный интеграл по достаточно простой области, но у меня более сложный случай и указанные формулы мне явно не подойдут.
Ну и так же я разбирался с методом Монте-Карло как таковым.

Для примера, мне надо взять такой интеграл
$\int\limits_{1}^{+\infty}dx \int\limits_{\arcsin(\sqrt{1/x})}^{\pi/2}\sin\theta d\theta \int\limits_{\arccos(\frac{1}{\sqrt{x}\sin\theta})}^{2\pi-\arccos(\frac{1}{\sqrt{x}\sin\theta})}d\varphi (-\cos\varphi\sin\theta+\sqrt{\frac{1}{x}})xe^{-\alpha x\cdot f(\theta,\varphi)}$

Мне бы хотелось применить метод Монте-Карло, потому что позже этот интеграл усложнится, он может стать 5-и кратным, и область интегрирования станет тоже сложнее. Можно ли как-то применить метод Монте-Карло для подобного интеграла?

PS, просьба модераторам, не могли бы вы переименовать тему в "Кратный интеграл по сложной области, как посчитать численно?" и перенести тему в раздел "численные методы". Заранее спасибо!

[AKM]

TheEnt в сообщении #511970 писал(а):

не могли бы вы ... перенести тему в раздел "численные методы"

Именные разделы --- не место для дискуссий*. Редкая тема (в процессе архивации) попадает в такой раздел.

*) Кроме нарушителей правил: они иногда там дискутируют.

[TheEnt]

AKM в сообщении #511989 писал(а):

Именные разделы --- не место для дискуссий*. Редкая тема (в процессе архивации) попадает в такой раздел.

Понял, спасибо!

[Sonic86]

TheEnt в сообщении #511970 писал(а):

Ещё раз, как звучит вопрос - как численно взять кратный интеграл со сложными пределами? Можно ли использовать метод Монте-Карло? Почти все, что я читал, применимо только для интегрирования по n-мерному параллелепипеду.

Можно через Монте-Карло (правда лично я не понимаю, на фига нужен Монте-Карло, если можно взять правильную решетку точек, на которых и считать в цикле значение функции).

TheEnt в сообщении #511970 писал(а):

Что я читал по моему вопросу - Демидович Б.П., Марон И.А. ... Калиткин Н.Н. - Численные методы

Есть еще Бахвалов. Довольно толстая книга, наверняка что-то найдется об этом. (Вот я тут, в частности, вижу параграф "Ускорение сходимости метода Монте-Карло". Да и вообще глава "Многомерные задачи")

TheEnt в сообщении #511970 писал(а):

Для примера, мне надо взять такой интеграл
$\int\limits_{1}^{+\infty}dx \int\limits_{\arcsin(\sqrt{1/x})}^{\pi/2}\sin\theta d\theta \int\limits_{\arccos(\frac{1}{\sqrt{x}\sin\theta})}^{2\pi-\arccos(\frac{1}{\sqrt{x}\sin\theta})}d\varphi (-\cos\varphi\sin\theta+\sqrt{\frac{1}{x}})xe^{-\alpha x\cdot f(\theta,\varphi)}$

А Вы предварительно переходили в сферические координаты? Если да, то Монте-Карло, конечно, можно, только переходить в сферические координаты незачем, если описание области будет проще - прога будет работать быстрее. Кроме того, будет необходимо избавиться от бесконечных пределов какой-нибудь специальной заменой (дробно-рациональной).

[_hum_]

Sonic86 в сообщении #512050 писал(а):

Можно через Монте-Карло (правда лично я не понимаю, на фига нужен Монте-Карло, если можно взять правильную решетку точек, на которых и считать в цикле значение функции).

Правильная решетка не оптимальная в плане скорости сходимости (слишком много в ней ненужных для получения искомого результата точек). При подходящем же выборе точек, например, из последовательности Холтона, говорят, можно получить скорость вплоть до $O(1/n^{1-\varepsilon}$) (а это даже быстрее, чем в обычном методе Монте-Карло, в котором скорость оценивается как $O(1/n^{1/2}$) ).

[Sonic86]

_hum_ в сообщении #512064 писал(а):

Правильная решетка не оптимальная в плане скорости сходимости (слишком много в ней ненужных для получения искомого результата точек). При подходящем же выборе точек, например, из последовательности Холтона, говорят, можно получить скорость вплоть до $O(1/n^{1-\varepsilon}$) (а это даже быстрее, чем в обычном методе Монте-Карло, в котором скорость оценивается как $O(1/n^{1/2}$) ).

Вооот оно что, спасибо А нам-то в универе говорили, вот тупо Монте-Карло и все. Даже прогу писали

Насчет неоптимальности правильной решетки по сравнению с простой генерацией случайных точек не верю (должно быть равноценно, но у Монте-Карло погрешность вероятностная, что хуже). Дадите ссылку на доказательство?

[_hum_]

Не совсем доказательство, но там, по-моему, есть ссылки на первоисточники:
И.М. Соболь. Численные методы Монте-Карло, М. 1973, глава 7 "Неслучайные точки в алгоритмах Монте-Карло"

[svv]

TheEnt, интересно записать Ваш интеграл в декартовых координатах. У меня тоже такое чувство, что Вы с них и начинали.
Что мне понятно: $x=r^2$. (Но далее $x$ везде в смысле декартовой координаты). Пределы по $\varphi$ выбраны так, чтобы выполнялось условие $x\leqslant 1$. Верхний предел $\theta=\pi/2$ дает условие $z\geqslant 0$.

А вот нижний предел по $\theta$ не очень понятен -- он означает, что $\rho=\sqrt{x^2+y^2} \geqslant 1$. Получается странное пересечение внешности единичного цилиндра и полупространства $x\leqslant 1$ (ну, и $z\geqslant 0$).

В математике, конечно, могут быть любые области, но у Вас ведь, скорее всего, физика...

Но вот если допустить, что Вы при наборе нечаянно написали $\arcsin(1/\sqrt{x})$ вместо $\arccos(1/\sqrt{x})$, тогда это давало бы гораздо более приятное условие $z\leqslant 1$.

[мат-ламер]

Можно предложить ещё идеи. 1) Двойной интеграл свести к последовательности одномерных интегралов. 2) Разбить область на треугольники и интерполировать интегрируемую функцию кубическим сплайном.

[TheEnt]

Все большое спасибо! Да, задача физическая, я свел её к взятию такого интеграла.

Sonic86 в сообщении #512050 писал(а):

А Вы предварительно переходили в сферические координаты? Если да, то Монте-Карло, конечно, можно, только переходить в сферические координаты незачем, если описание области будет проще - прога будет работать быстрее. Кроме того, будет необходимо избавиться от бесконечных пределов какой-нибудь специальной заменой (дробно-рациональной).

Интеграл конечно можно записать и в декартовых координатах, но дело в том, что я написал здесь его немного в упрощенном виде для примера, подынтегральное выражение ещё нужно будет дополнительно умножить на функцию, зависящую от всех трех переменных, в декартовых переменных будет неудобно. Да, бесконечный предел надо убрать, спасибо, что подсказали!

svv в сообщении #512125 писал(а):

TheEnt, интересно записать Ваш интеграл в декартовых координатах. У меня тоже такое чувство, что Вы с них и начинали.
Что мне понятно: . (Но далее везде в смысле декартовой координаты). Пределы по выбраны так, чтобы выполнялось условие . Верхний предел дает условие .

А вот нижний предел по не очень понятен -- он означает, что . Получается странное пересечение внешности единичного цилиндра и полупространства (ну, и ).

В математике, конечно, могут быть любые области, но у Вас ведь, скорее всего, физика...

Но вот если допустить, что Вы при наборе нечаянно написали вместо , тогда это давало бы гораздо более приятное условие .

Область именно такая, я записал все верно. Она определяется условием $(\cos\varphi\sin\theta-\sqrt{\frac{1}{x}})<0$ (вроде знак не перепутал). Можно интерпретировать так - проекция единичного вектора на ось $Ox$ должна быть меньше, чем $\sqrt{\frac{1}{x}}$, и в крайних положениях по $x$ получается интегрирование либо по всей сфере, либо по полусфере.

Вообще конкретный вид интеграла сейчас не важен, скорее всего он будет усложнен. Тут происходит интегрирование по координатам импульсного пространства, я собираюсь ввести ещё зависимость от обычных пространственных координат. Поэтому я хочу делать все методом Монте-Карло, делать сетку для 5-6 кратного интеграла не хочется, да и сходимость будет хуже (если я не ошибаюсь...)

Если можно делать методом Монте-Карло, как я должен выбрасывать случайные точки? Должен ли я заключить свою область в n-мерный параллелепипед, распространить подынтегральную функцию на новую область с помощью функции Хевисайда и выбрасывать точки в нем, или я должен выбрасывать точки только в области интегрирования?

Я сегодня попробовать сделать элементарную задачу, я взял интеграл $\int\limits_{0}^{1}rdr\int\limits_{0}^{r^2}dz\int\limits_{0}^{2\pi}d\varphi=\frac{\pi}{2}$ методом Монте-Карло, сведя его к $2\pi\int\limits_{-1}^{1}dr\int\limits_{-1}^{1}dz\int\limits_{-1}^{1}d\varphi r\theta(r)\theta(z)\theta(r^2-z)\theta(\varphi)$, и с увеличением числа точек точность увеличивалась. Когда я брал этот интеграл по кубатурной формуле Симпсона, при увеличении количества точек точность не увеличивалась.