Несобственные кратные интегралы в Maple

[Mikhail_K]

Здравствуйте.
Такая проблема: нужно научиться считать в Maple несобственные кратные интегралы, численно.
Например, я пишу

int(int(sin(x*y+2)/sqrt(x^2+y^2),x=-1..1),y=-1..1);

Просто выдаёт символ интеграла, не считает.
Пробую по-другому:

int(int(sin(x*y+2)/sqrt(x^2+y^2),x=-1..1),y=-1..1,numeric);

numeric - чтобы считал численно. Выдаёт очень странный ответ про рекурсию:

Error, (in tools/map) too many levels of recursion

Откуда здесь рекурсия-то?
Пробовал писать numeric у обоих интегралов, всё равно не помогает.
Видел на форуме совет писать так:

MyInt:=Int(Int(sin(x*y+2)/sqrt(x^2+y^2),x=-1..1),y=-1..1);
evalf[8](MyInt);

То же самое.
Ещё слышал, что для несобственных надо применять опцию continuous. Тоже не считает.
Помогите, что делать?

[GAA]

Если точность не важна, то попробуйте так

> evalf(Int(sin(x*y+2)/sqrt(x^2+y^2), [x=-1..1, y=-1..1], method=_MonteCarlo, epsilon=0.00005));
                                   6.222174505

Если точность важна, то попробуйте перебрать имеющиеся методы. Если не получится (напишите, какие перебрали методы), то придется разбираться в математике и программировать.

[Otta]

Mikhail_K в сообщении #848410 писал(а):

int(int(sin(x*y+2)/sqrt(x^2+y^2),x=-1..1),y=-1..1);

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

[Mikhail_K]

Otta в сообщении #848427 писал(а):

По существу проблемы: ее источник совершенно понятен, Вы считаете повторный интеграл вместо кратного, а он не существует. Уже внутренний расходится.

Вы ошибаетесь. Внутренний интеграл расходится только при y=0. То есть, внешний интеграл берётся от функции, которая обращается в бесконечность при y=0. Другими словами, внешний интеграл просто несобственный, вот и всё.

[Otta]

Внутренний тоже несобственный. Когда $y=0$. И при этом значении $y$ он превращается в $\int_{-1}^1\frac{\sin 2}{|x|}\,dx$. Сходится?

[Mikhail_K]

Внутренний РАСходится при y=0, и только при y=0.
Это значит, что внутренний интеграл есть функция от y, которая обращается в бесконечность при y=0, и принимает конечные значения при других y. Такую функцию можно интегрировать по y. Таким образом, повторный интеграл (именно повторный) может существовать, несмотря на расходимость внутреннего при y=0.

[Otta]

Mikhail_K в сообщении #848437 писал(а):

Таким образом, повторный интеграл (именно повторный) может существовать, несмотря на расходимость внутреннего при y=0.

Это Вы знаете, это я знаю. Вопрос: знает ли это машина? Иными словами, не получается ли, что утыкаясь в неприличное поведение интеграла при одном значении $y$, она просто не знает, что с этим делать? Какой заложен алгоритм вычисления?

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

[Mikhail_K]

Синус наверху я поставил затем, что без него этот интеграл берётся, и Maple его просто берёт, а не считает численно.

[Otta]

А от $\frac{1}{(x^4+y^4)^{1/3}}$ что получается?

[GAA]

Otta
Maple 7, Classic Worksheet Maple 12, Classic Worksheet Maple 15

> evalf(Int(Int(1/(x^2+y^2)^(1/2), x=-1..1), y=-1..1));
                        Вычислить не может
> evalf(Int(Int(1/(x^4+y^4)^(1/3), x=-1..1), y=-1..1));
                        Float(infinity)

Classic Worksheet Maple 12

> evalf(Int(1/(x^4+y^4)^(1/3), [x=-1..1, y=-1..1], method=_MonteCarlo, epsilon=0.0005)); 
                                  11.37465800
> evalf(Int(1/(x^2+y^2)^(1/2), [x=-1..1, y=-1..1], method=_MonteCarlo, epsilon=0.0005));
                                    7.050293210

Перейдя в полярную систему координат, в случае квадратного корня в Maple 12 получим

 > 2*(evalf(Int(Int(1, rho = 0..1/cos(phi)), phi = -Pi/4..Pi/4)) + evalf(Int(Int(1, rho=0..1/sin(phi)), phi = Pi/4..3*Pi/4)));
                                     7.050988696

А в случае кубического корня в Maple 12

 > f:= 1/(cos(phi)^4+sin(phi)^4)^(1/3): g:= 1/rho^(1/3):
> 2*(evalf(Int(f*Int(g, rho = 0..1/cos(phi)), phi = -Pi/4..Pi/4)) + evalf(Int(f*Int(g, rho=0..1/sin(phi)), phi = Pi/4..3*Pi/4)));
                                       11.39721950 

Аналогично, в случае исходного интеграла

 > f:= sin(rho^2*sin(phi)*cos(phi)+2):
> 2*(evalf(Int(Int(f, rho = 0..1/cos(phi)), phi = -Pi/4..Pi/4)) + evalf(Int(Int(f, rho=0..1/sin(phi)), phi = Pi/4..3*Pi/4)));
                                       6.222736680

Результат близок к значению, полученному при помощи метода _MonteCarlo.
Попробую грубо оценить относительную погрешность значения, полученного во втором сообщении темы методом _MonteCarlo. Для этого увеличу точность вычислений и буду вычисленное с такой точностью значение считать «точным».

> Digits:=20;f:= sin(rho^2*sin(phi)*cos(phi)+2):
> ff:=2*(evalf(Int (Int(f, rho = 0..1/cos(phi)), phi = -Pi/4..Pi/4)) + evalf(Int(Int(f, rho=0..1/sin(phi)), phi = Pi/4..3*Pi/4))):
> eps:= abs((ff - 6.222174505)/ff);
              eps := .000090342106339135968349

Увы, возникает подозрение, что метод _MonteCarlo не обеспечил заявленную погрешность вычислений. Хоть и не сильно.

Mikhail_K
Фичи и баги Maple — это необозримая тема. Я так и не понял, что планируется обсуждать в данной ветке.

-- Sat 12.04.2014 07:46:44 --

Документация Online по numeric integration Maple 17.

-- Соединены подряд идущие сообщения c перерывом три с половиной часа --

GAA в [url=http://dxdy.ru/post848561.html#p848561]сообщении #848561 писал(а):

Maple 7, Classic Worksheet Maple 12, Classic Worksheet Maple 15

> evalf(Int(Int(1/(x^2+y^2)^(1/2), x=-1..1), y=-1..1));
                        Вычислить не может
> evalf(Int(Int(1/(x^4+y^4)^(1/3), x=-1..1), y=-1..1));
                        Float(infinity)

Это я по невнимательности не в ту программу данные вел [на сервере маленький gray scale монитор]. На самом деле, maple 12 и maple 15 пишут, что не могут вычислить evalf(Int(Int(1/(x^4+y^4)^(1/3), x=-1..1), y=-1..1));. Т.е. эту дырку Maple 7 как будто починили. Остальное перепроверил, подтверждается.

[Mikhail_K]

На самом деле, как оказалось, несобственные интегралы в Maple можно считать очень просто.
Именно, вместо

int(sin(x*y+2)/sqrt(x^2+y^2),[x=-1..1,y=-1..1],numeric);

писать

eps:=10^(-12);
int(sin(x*y+2)/sqrt(x^2+y^2),[x=-1..-eps,y=-1..1],numeric)+
int(sin(x*y+2)/sqrt(x^2+y^2),[x=eps..1,y=-1..1],numeric);

Точность замечательная

[Otta]

Mikhail_K в сообщении #848664 писал(а):

int(sin(x*y+2)/sqrt(x^2+y^2),[x=-1..1,y=-1..1],numeric);

Интересно, а эта команда что, не работает? Ибо второе предложение - некий изврат, основанный на использовании определения.

[Mikhail_K]

Otta в сообщении #848683 писал(а):

Интересно, а эта команда что, не работает? Ибо второе предложение - некий изврат, основанный на использовании определения.

Не работает. Долго думает, затем выдаёт ошибку.
Если подынтегральная функция непрерывна, такая команда работает, считает мгновенно. А вот уходящие на бесконечность функции Maple, видать, не любит.
Одномерные несобственные интегралы считает, а вот такие кратные несобственные - нет.
Только не соглашусь, что это изврат. Maple умеет считать собственные интегралы. Вычисление несобственного интеграла сводится к вычислению собственных по определению.

А с методом Монте-Карло перестаёт работать при требовании большой точности. Здесь же - какая угодно точность. И скорость великолепная.

[GAA]

На всякий случай добавлю: вызов
int(sin(x*y+2)/sqrt(x^2+y^2),[x=-1..1,y=-1..1],numeric);
эквивалентен (по справке) вызову
int(int(sin(x*y+2)/sqrt(x^2+y^2),x=-1..1),y=-1..1);
(Более того, раньше первого варианта записи двойного интеграла не было и следовало двойной интеграл задавать только в виде повторного. Так было еще в Maple 7. Точно с какой версии появился первый тип записи я не скажу. Но в 12 версии он точно есть.)

Об "изврате". Стандартный подход: заменяем несобственный интеграл собственным + малая добавка. В данном случае выбрасываем полосу, в которой находится особая точка. Разве не этому учат студентов в курсе «Высшей математики» или «Численных методов»? Недостатком этого подхода является необходимость оценить ширину полосы, которую выбрасываем. И с увеличением точности вычисления интеграла eps следует пересчитывать.
Преимуществом перехода в полярную систему координат будет большая скорость вычислений. Однако, если нам точность не важна и например 10 знаков достаточно, то эти тонкости не имеют значения.

Maple некоторые кратные несобственные интегралы вычисляет, некоторые нет. Заранее не угадаешь. Каждый раз приходиться разбираться индивидуально.

Редактирование: numeric лишний, зачеркнул.

[Otta]

Так может, и считать его сразу по определению? и по обеим переменным выкалывать окрестность?
Кстати, Вы не ответили, что получается для функции

Otta в сообщении #848445 писал(а):

$\frac{1}{(x^4+y^4)^{1/3}}$

. Но это уже выяснили, а "по определению" что получается?

А, вот,

GAA в сообщении #848692 писал(а):

На всякий случай добавлю: вызов
int(sin(x*y+2)/sqrt(x^2+y^2),[x=-1..1,y=-1..1],numeric);
эквивалентен (по справке) вызову
int(int(sin(x*y+2)/sqrt(x^2+y^2),x=-1..1),y=-1..1);

я так и подозревала. :(