L'intégration d'une fonction avec singularites en utilisant la routine de quad scipy

Question

J'utilise la fonction quad de scipy.integrate v0.19.1 d'intégrer des fonctions avec une racine carrée comme singularité à chaque extrémité de l'intervalle d'intégration comme par exemple

In [1]: quad(lambda x: 1/sqrt(1-x**2), -1, 1) 

(I utiliser la fonction sqrt de numpy v1.12.0) qui donne immédiatement le pi résultat correct:

Out[1]: (3.141592653589591, 6.200897573194197e-10) 

Selon la documentation du quad fonction le mot-clé points doit être utilisé pour indiquer l'emplacement des singularités ou discontinuités du integrand, mais si j'indiquer les points [1, -1] où le integrand ci-dessus est singluar je reçois un avertissement et nan comme résultat:

In [2]: quad(lambda x: 1/sqrt(1-x**2), -1, 1, points=[-1, 1]) 

RuntimeWarning: divide by zero encountered in double_scalars 
IntegrationWarning: Extremely bad integrand behavior occurs at some 
points of the integration interval. 

Out[2]: (nan, nan) 

Can quelqu'un clarifier, pourquoi quad produit ces problèmes si les singularités de l'intégrande sont spécifiées et fonctionne bien si ces points ne sont pas indiqués?

EDIT: Je pense que j'ai trouvé la bonne solution pour ce problème. Pour le cas où quelqu'un rencontre les autres problèmes similaires que je veux rapidement partager mes conclusions:

Je vous souhaite intégrer une fonction de la forme f(x)*g(x) avec une fonction lisse f(x) et g(x) = (x-a)**alpha * (b-x)**beta, où a et b sont les limites d'intégration et g(x) a des singularités à ces limites si alpha, beta < 0, alors vous êtes supposé simplement intégrer f(x) en utilisant g(x) comme fonction de pondération. Pour la routine quad, cela est possible en utilisant les arguments weight et wvar. Avec ces arguments, vous êtes également capable de gérer des singularités de différents types et des comportements oscillatoires problématiques. La fonction de pondération g(x) définie ci-dessus peut être utilisée en définissant weight='alg' et en utilisant wvar=(alpha, beta) pour spécifier les exposants dans g(x).

Depuis 1/sqrt(1-x**2) = (x+1)**(-1/2) * (1-x)**(-1/2) je peux maintenant gérer l'intégrale comme suit:

In [1]: quad(lambda x: 1, -1, 1, weight='alg', wvar=(-1/2, -1/2)) 
Out[1]: (3.1415926535897927, 9.860180640534107e-14) 

qui donne la bonne réponse pi avec une très grande précision, peu importe si j'utilise l'argument points=(-1, 1) ou non (qui, dans la mesure Je comprends maintenant, ne devrait être utilisé, si les singularités/discontinuités ne peuvent pas être traitées en choisissant une fonction de pondération appropriée).

Answer 1

Le paramètre points est destiné aux singularités/discontinuités qui se produisent au sein de l'intervalle d'intégration. Il n'est pas destiné aux extrémités de l'intervalle. Ainsi, dans votre exemple, la version sans points est la bonne approche. (Je ne peux pas identifier ce qui ne va pas quand les points d'extrémité sont inclus dans points sans plonger dans le code FORTRAN que SciPy enveloppe.)

Comparer avec l'exemple suivant, où se produisent dans l'singularites intervalle d'intégration:

>>> quad(lambda x: 1/sqrt(abs(1-x**2)), -2, 2) 
(inf, inf) 
>>> quad(lambda x: 1/sqrt(abs(1-x**2)), -2, 2, points = [-1, 1]) 
(5.775508447436837, 7.264979728915932e-10) 

Ici l'inclusion de points est appropriée, et donne le bon résultat, alors que sans points la sortie est sans valeur.