могу ли я доверять этому числовому тройному интегралу из Matlab?

15

Вычислительные люди науки:

Первоначально я разместил этот вопрос на Math Stack Exchange, и кто-то сказал, что я мог бы получить «гораздо лучшие» ответы здесь:

Я новичок в численных методах и Matlab. Я пытаюсь оценить следующую сумму двух тройных интегралов (она, очевидно, может быть записана более просто, но вы все еще не можете оценить ее символически (?)). У меня возникают проблемы с тем, чтобы работал здесь, поэтому я неохотно разбил его на куски: я хочу найти суммуLATЕИкс

2((1/0,3)-1)2(11/0,31р10р1-р0F1(р0,р1,T)ехр(-(0,3)2T24)dTdр0dр1),

и

2((1/0,3)-1)2(11/0,31р1р1-р0р1+р0F2(р0,р1,T)ехр(-(0,3)2T24)dTdр0dр1),

где

F1(р0,р1,T)знак равноT2р03*(0,3)32р13π

и

F2(р0,р1,T)знак равно(0,3)3π3/2(р0+р1-T)4(T2+2T(р0+р1)-3(р1-р0)2)2288(43πр03)(43πр13),

РЕДАКТИРОВАТЬ (2 марта 2013 г.): Кто-то ответил, что они заставили Mathematica сделать интегралы символически. Я только попытался сделать это (с упрощенными версиями интегралов), и Mathematica могла сделать только два внешних из первого, и остановилась на втором. Буду признателен за помощь. Вот что я сделал.

Я пытался оценить

121р20р2-р1р13T2ехр(-T2)р23dTdр1dр2
через

Интегрировать [r1 ^ 3 / r2 ^ 3 * t ^ 2 * Exp (-t ^ 2), {t, 0, r2 - r1}, {r1, 1, r2}, {r2, 1, 2}]

и Mathematica возвращается (у меня были проблемы с здесь, потому что результат длинный. Я разбил его на два уравнения. Если кто-нибудь знает хороший способ показать это, пожалуйста, сообщите мне):LATЕИкс

12164р22е-1-р22(2е2р2(25+р2(19+2р2(1+р2)))-

е1+р22(32р2(2+р22))+π(11+4р22(9+р22))приусадебный участок[1-р2])dр2.

Тогда я попытался оценить

121р2р2-р1р2+р1...

...ехр(-T2)(р1+р2-T)4(T2+2T(р1+р2)-3(р2-р1)2)2р13р23dTdрdр2

с помощью

Интегрировать [(r1 + r2 - t) ^ 4 * (t ^ 2 + 2 * t * (r1 + r2) - 3 * (r2 - r1) ^ 2) ^ 2 * Exp [-t ^ 2] / r1 ^ 3 / r2 ^ 3, {r2, 1, 2}, {r1, 1, r2}, {t, r2-r1, r2 + r1}]

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

[Конец редактирования 2 марта]

Я использовал команду "triplequad" от Matlab без дополнительных опций. Я управлял переменными пределами интеграции с помощью тяжелых функций, потому что я не знал другого способа сделать это. Матлаб дал мне . +0,007164820144202

Я знаю, что Matlab - хорошее программное обеспечение, но я слышал, что числовые тройные интегралы трудно сделать точно, и математики должны быть настроены скептически, поэтому я хочу как-то проверить точность этого ответа. Интегралы дают ожидаемую ценность определенного эксперимента (если кто-то хочет, я могу отредактировать этот вопрос, чтобы описать эксперимент): я реализовал эксперимент в Matlab, используя соответственно случайно сгенерированные числа, миллион раз, и усреднил результаты. Я повторил этот процесс четыре раза. Вот результаты (я прошу прощения, если я неправильно использовал слово «испытание»):

Пробная версия 1:0.007133292603256

Пробная версия 2:0.007120455071989

Пробная версия 3:0.007062595022049

Пробная версия 4:0.007154940168452

Пробная версия 5:0.007215000289130

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

Так может ли кто-нибудь сказать мне, могу ли я доверять результату «тройной четверки» здесь, и при каких обстоятельствах можно доверять ему вообще?

Одним из предложений, которое я получил на Math Stack Exchange, было попробовать другое программное обеспечение, такое как Mathematica, Octave, Maple и SciPy. Это хороший совет? Люди действительно выполняют числовую работу в Mathematica и Maple? Octave - это своего рода клон Matlab, поэтому можно ли предположить, что он использует те же алгоритмы интеграции? Я даже не слышал о SciPy раньше и был бы признателен за любые мнения по этому поводу.


ОБНОВЛЕНИЕ: Кто-то из Math Stack Exchange сделал это в Maple и получил . Это согласие с тремя значительными цифрами. Это хороший знак.0.007163085468

Кроме того, я был бы признателен за предложения о том, как ввести длинное многострочное выражение в в Stack Exchange. Можете ли вы использовать "выровненную" среду здесь? Я пытался, и я не мог заставить его работать.LATEX

Стефан Смит
источник
2
Ваши результаты моделирования прекрасно согласуются с численным значением , возвращаемым Matlab: их среднее только стандартные ошибки меньше , чем возвращается Matlab. FWIW, Mathematica возвращает . Он также может символически оценивать эти интегралы в терминах полиномов и функций ошибок. 0.007137261.110.00716308537
whuber
@whuber Спасибо. Я мог поклясться, что попробовал это символически в Maple, а Maple не смог. Я попробую еще раз в Maple, и если это не сработает, я попробую это в Mathematica. Кстати, я сделал аналогичный интеграл в Maple, и я получил огромный символический ответ. Казалось, что это сумма и разница очень больших чисел, общая сумма которых была довольно мала. Я подозреваю, что ошибка округления была, вероятно, в окончательном ответе. В такой задаче вы должны использовать символический ответ или просто сделать числовой интеграл?
Стефан Смит
Символьные ответы имеют то преимущество, что они представляют собой комбинации функций, которые (часто) могут эффективно вычисляться с произвольной точностью. Как правило, символическое решение также позволяет быстро пересчитывать при изменении параметров. По этим причинам часто стоит искать символическое решение.
2011 г.,
@whuber: Я попытался сделать некоторые по существу эквивалентные интегралы (изменив некоторые константы и удалив некоторые мультипликативные константы) в Mathematica, и Mathematica могла выполнять только два внешних интегрирования первого интеграла, и, кажется, остановилась на втором. Я разместил свой код и результаты выше.
Стефан Смит
1
Отредактируйте 2 марта: уменьшив тройной интеграл символически до единого (в первой половине ваших интегралов), вы добились многого. Интегрант очень хорошо ведет себя и может быть численно интегрирован с чрезвычайно высокой точностью за доли секунды.
2013 г.

Ответы:

9

108

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

Если вы захотите сделать это с помощью SciPy, вам также потребуется прибегнуть к вложенной адаптивной квадратуре Гаусса, используя базовые подпрограммы Quadpack (Piessens et al.). В Октаве у вас будет такой же подход. И я бы не слишком удивлен , если Matlab также использует Quadpack в качестве квадратурного двигателя (так как это ссылка).

GertVdE
источник
@GretVdE: Спасибо за информацию. Сначала я попытался символически оценить интеграл, и Maple не смог этого сделать (поэтому, возможно, это невозможно, используя стандартные функции), поэтому я попросил Maple сделать это численно. Я не знаю, какой алгоритм он использовал.
Стефан Смит
@StefanSmith: вы можете узнать, установив InfoLevel в Maple: infolevel[`evalf/int`] := 4. Вы уверены, что Mape не может найти закрытое решение? Интеграл не кажется слишком сложным. Не могли бы вы опубликовать свой кленовый лист где-нибудь?
GertVdE
@StefanSmith: я бы опубликовал код Maple в вопросе выше.
GertVdE
Я не могу заставить Maple работать в моей системе прямо сейчас, но я попытался использовать эквивалентные интегралы в Mathematica, и Mathematica выполнила только два внутренних из первого тройного интеграла и остановилась на втором тройном интеграле. Пожалуйста, смотрите отредактированный вопрос.
Стефан Смит