Czy ktoś może wyjaśnić następujące zachowanie, które obserwuję podczas integracji CHEBYSHEVE Funkcja wagi przy użyciu 3 różnych procedur i 2 różne reprezentacje wykładnika? Oczekiwana odpowiedź brzmi PI w każdym przypadku:

from scipy.integrate import quadrature, quad, fixed_quad

print fixed_quad(lambda x: 1/(1 - x**2)**(1/2), -1, 1)
print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

print quadrature(lambda x: 1/(1 - x**2)**(1/2), -1, 1) 
print quadrature(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

print quad(lambda x: 1/(1 - x**2)**(1/2), -1, 1)
print quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1)

Co daje

(2.0000000000000009, None)
(2.8254100794589787, None)

(1.9999999999999996, 4.4408920985006262e-16)
/usr/lib/python2.7/dist-packages/scipy/integrate/quadrature.py:168: AccuracyWarning:    maxiter (50) exceeded. Latest difference = 6.965869e-04
 AccuracyWarning)
(3.107110439388189, 0.00069658693569163432)

(2.0, 2.220446049250313e-14)
(3.141592653589564, 6.200200353134733e-10)

Po pierwsze, jak widać, odpowiedź zależy od tego, czy 1/2 lub 0,5 jest podana w wykładniku: Dlaczego tak jest?

Po drugie, wynik zależy od wybranej procedury; Czy ktoś może wyjaśnić, dlaczego QuadPack Fortrana dostaje właściwą odpowiedź, ale czworokraturze Gaussiańskie przegapią wynik całkowicie?

Notatka : Podobne pytanie jest adresowane Oto, ale nie dotyczy dwóch określonych pytań powyżej.

0
MaviPranav 13 sierpień 2014, 19:15

2 odpowiedzi

Najlepsza odpowiedź

Ponieważ @Gjdanis wskazuje, w Pythonie 2.7, {X0}} 0 (chyba że zawiera from __future__ import division w kodzie).

Twoja integranda ma osobliwości na 1 i -1. fixed_quad i quadrature Wykonaj kwadraturę Gaussa z funkcją wagi w(x) = 1, więc te osobliwości nie są dobrze obsługiwane.

fixed_quad nie jest adaptacyjny (stąd nazwa). Zamówienie domyślne to 5. Musisz zwiększyć zamówienie (przez wiele), aby uzyskać rozsądne przybliżenie:

In [179]: print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1, n=100)
(3.124265558250825, None)

In [180]: print fixed_quad(lambda x: 1/(1 - x**2)**(0.5), -1, 1, n=2000)
(3.1407221810853478, None)

quadrature Wystarczy wzywać fixed_quad z rosnącym porządkiem (do maksymalnego argumentu {X2}}), aż różnica między kolejnymi integralnymi szacunkami jest wystarczająco mała. Ostrzeżenie, które drukowane mówi, że osiągnięto maksymalne zamówienie bez spełnienia pożądanej tolerancji błędów. Domyślne dla maxiter wynosi 50; Musisz zwiększyć maxiter, aby uzyskać lepsze wyniki. Na przykład tutaj jest wynik z maxiter=200:

In [2]: print quadrature(lambda x: 1/(1 - x**2)**(0.5), -1, 1, maxiter=200)
/Users/warren/local_scipy/lib/python2.7/site-packages/scipy/integrate/quadrature.py:183: AccuracyWarning: maxiter (200) exceeded. Latest difference = 4.353464e-05
  AccuracyWarning)
(3.1329074742380407, 4.3534643496823122e-05)

Jeśli użyjesz maxiter, powinieneś także uniwać sędziowe użycie miniter. quadrature naiwnie zaczyna się w miniter i zwiększa zamówienie o 1, aż oszacowanie błędów jest wystarczająco małe lub maxiter zostanie osiągnięte.

Aby dowiedzieć się więcej o tym, jak fixed_quad i quadrature działa, spójrz na kod źródłowy: https://github.com/scipy/scipy/blob/master/scipy/integrate/quadrature.py

Jak wskazujesz, quad jest {} {{{x0} Ten kod jest znacznie bardziej wyrafinowany niż prosta czworokatura Gaussa używana przez fixed_quad i quadrature.

1
Warren Weckesser 13 sierpień 2014, 18:20

W Pythonie 2.7 (który używasz) Dział Integer jest domyślny. Oznacza to, że 1/2 oceni się do 0. Jeśli chcesz użyć podziału pływającego punktu, dodaj from __future__ import division na górę kodu.

3
rookie 13 sierpień 2014, 15:17