2009-02-24 10 views
4

Я пытаюсь численно интегрировать произвольную (известную, когда я код) функцию в моей программе с использованием численных методов интеграции. Я использую Python 2.5.2 вместе с пакетом численных интеграций SciPy. Для того, чтобы получить чувство для этого, я решил попытаться интегрировать грех (х) и наблюдал этого behavior-Может кто-нибудь объяснить, почему scipy.integrate.quad дает разные результаты для одинаково длинного диапазона при интеграции sin (X)?

>>> from math import pi 
>>> from scipy.integrate import quad 
>>> from math import sin 
>>> def integrand(x): 
...  return sin(x) 
... 
>>> quad(integrand, -pi, pi) 
(0.0, 4.3998892617846002e-14) 
>>> quad(integrand, 0, 2*pi) 
(2.2579473462709165e-16, 4.3998892617846002e-14) 

Я считаю такое поведение странного, потому что -
1. В обычной интеграции, интегрируя по полному циклу дает нуль.
2. При численном интегрировании это (1) не обязательно имеет место, поскольку вы можете просто быть , приближаясь к общей площади под кривой.

В любом случае, если предположить, что 1 имеет значение True или предполагается, что значение 2 истинно, я считаю поведение непоследовательным. Либо оба интегрирования (от -pi до pi и от 0 до 2 * pi) должны возвращать 0,0 (первое значение в кортеже является результатом, а второе - ошибкой) или return 2.257 ...

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

В любом случае, в моем конечном приложении, я планирую использовать вышеуказанный метод для нахождения длины дуги функции. Если у кого-то есть опыт в этой области, пожалуйста, сообщите мне о лучшей политике для этого в Python.

Edit
Примечание
У меня уже есть первые дифференциальные значения во всех точках в диапазоне, хранящихся в массиве.
Текущая ошибка допускается.
Конец примечания

Я прочитал Wikipaedia на этом. Как отметил Димитрий, я буду интегрировать sqrt (1 + diff (f (x), x)^2), чтобы получить длину дуги. Я хотел спросить, есть ли лучшее приближение/Лучшая практика (?)/Более быстрый способ сделать это. Если требуется больше контекста, я отправлю его отдельно/пост контекста, как вы пожелаете.

+0

Вот сообщение в блоге http://algorithmist.wordpress.com/2009/01/05/quadratic-bezier-arc-length/, которое может показаться вам полезным. Ссылка из блога http://www.algorithmist.net/crlength.html – jfs

+0

Спасибо. Я посмотрю ссылки. – batbrat

ответ

9

Функция quad является функцией из старой библиотеки Fortran. Он работает, судя по плоскостности и наклону функции, которую он интегрирует, как обрабатывать размер шага, который он использует для численного интегрирования, чтобы максимизировать эффективность. Это означает, что вы можете получить несколько разные ответы от одного региона к другому, даже если они аналитически одинаковы.

Без сомнения, обе интеграции должны возвращать ноль. Возвращение того, что составляет 1/(10 триллионов), довольно близко к нулю! Небольшие различия обусловлены тем, как quad перекатывается за sin и меняет размеры шага. Для вашей запланированной задачи вам понадобится quad.

EDIT: За что вы делаете Я думаю, quad в порядке. Это быстро и довольно точно. Мое окончательное утверждение использует его с уверенностью, если вы не найдете то, что действительно пошло наперекосяк. Если он не возвращает бессмысленный ответ, то, вероятно, он работает нормально. Не беспокойся.

+0

+1. Спасибо, что объяснили о квадранте, а также отметили, что квадрат достаточен. – batbrat

4

Этот вывод кажется правильным для меня, так как здесь у вас есть абсолютная оценка погрешности. Интегральное значение sin (x) действительно должно иметь значение нуля для полного периода (любой интервал длины 2 * pi) как в обычном, так и в цифровом интегрировании, и ваши результаты близки к этому значению.
Чтобы оценить длину дуги, вы должны вычислить интеграл для функции sqrt (1 + diff (f (x), x)^2), где diff (f (x), x) является производной от f (x). См. Также Arc length

+0

Спасибо за ответ. Я знаю о вычислении длины дуги. Я также прочитал статью wikipaedia. Я отредактирую свой пост, чтобы это отразить. – batbrat

+0

Еще раз спасибо! Я не заметил, что ответ был ниже предела ошибки. – batbrat

6

Я думаю, что это, вероятно, машинная точность, поскольку оба ответа фактически равны нулю.

Если вы хотите получить ответ от первоисточника, я бы разместить этот вопрос на scipy discussion board

+0

Спасибо за ссылку на дискуссионный форум. – batbrat

6

Я бы сказал, что число O (10^-14) эффективно нулю. Какова ваша терпимость?

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

Это может быть только природа чисел с плавающей точкой: "What Every Computer Scientist Should Know About Floating Point Arithmetic".

+0

Я могу терпеть ошибки более высокой величины, чем 10^-14 – batbrat

3
0.0 == 2.3e-16 (absolute error tolerance 4.4e-14) 

Оба ответа являются одинаковыми и правильно т.е., нулевой в пределах заданного допуска.

2

Отличие состоит в том, что sin (x) = - sin (-x) точно даже в конечной точности. В то время как конечная точность дает только sin (x) ~ sin (x + 2 * pi). Конечно, было бы неплохо, если бы квадрат был достаточно умен, чтобы понять это, но на самом деле не имеет никакого способа узнать априори, что интеграл в течение двух интервалов, которые вы даете, эквивалентен или что первый - лучший результат.

 Смежные вопросы

  • Нет связанных вопросов^_^