0

I am using quadpy to integrate a function in python.

Function

import numpy as np
T = 2*np.pi 

def ex1(t):
    return np.where(np.logical_and((t%T>=0), (t%T<np.pi)), t%T, np.pi)

The function is periodic, this is its plot:

x = np.linspace(0, 6*T, 1000)
plt.plot(x, ex1(x))
plt.grid(True)
plt.show()

enter image description here

Problem

I am trying to integrate this function:

from scipy.integrate import quad
import quadpy

print(quadpy.quad(ex1, 0, 3))
print(quad(ex1, 0, 3))

produces

(array(4.5), array(1.41692995e-19))
(4.5, 4.9960036108132044e-14)

On an interval from 0 to 3, everything works fine. However If I increae the interval to e.g. 4, scipy still works.

print(quad(ex1, 0, 4))

produces

(7.631568411183528, 1.0717732083155035e-08)

but

print(quadpy.quad(ex1, 0, 4))

produces

IntegrationError: Tolerances (abs: 1.49e-08, rel: 1.49e-08) could not be reached with the given max_num_subintervals (= 50).

Questions

  • How do I prevent this error? I tried adding an argument called max_num_subintervals but that did not seem to work.
  • Am I using quadpy correctly for what I am trying to do? I have begun using it since I wanted to take derivatives of complex numbers which scipy does not support, and I would like to have a one-size-fits-all solution, thus using quadpy for these easier examples that scipy would be enough.
charelf
  • 3,103
  • 4
  • 29
  • 51
  • Sounds like a bug that should be reported on github/quadpy. https://github.com/nschloe/quadpy/issues – Nico Schlömer Apr 14 '20 at 08:11
  • Hi, thanks for noticing my problem, I was hoping you would see my question! I have posted an answer that kind of prevents my problem, however I am unsure whether it is correct and why such a high tolerance is needed when scipy does not. Maybe you know more. – charelf Apr 14 '20 at 08:17
  • Like I said: Bug report! (This way I always see it.) – Nico Schlömer Apr 14 '20 at 08:17
  • Ok, I will post it. edit: Done https://github.com/nschloe/quadpy/issues/255 – charelf Apr 14 '20 at 08:18

2 Answers2

1

This was a quadpy bug after all. Fixed now:

import numpy as np
from scipy.integrate import quad
import quadpy

T = 2 * np.pi


def ex1(t):
    return np.where(np.logical_and((t % T >= 0), (t % T < np.pi)), t % T, np.pi)


print(quadpy.__version__)
print(quadpy.quad(ex1, 0, 4))
print(quad(ex1, 0, 4))

0.14.7
(array(7.63156841), array(1.13827355e-16))
(7.631568411183528, 1.0717732083161812e-08)

You're probably using an outdated quadpy version.

Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
  • I get the exact same results executing this code. Please try with `print(quadpy.quad(ex1, 0, 4))`, mine works with an interval from 0 to 3, but not 0 to 4. Other than that, quadpy version I have is the same. – charelf Apr 14 '20 at 18:12
  • Thanks for the update to quadpy 14.7, which solves the issue for intervals from e.g.0 to 4. However please try now with an interval from 0 to 7, or please see my comment in the github issue. Sorry for the inconveniences I'm causing, and thanks for your work! – charelf Apr 15 '20 at 08:27
0

After some more experimenting I have found that quadpy's quad method takes the same argument as scipys quad method, which can be found here:

Using the epsabs, epsrel, limit optional arguments I can thus prevent the error:

print(quadpy.quad(ex1, 0, 4, epsabs=1e-1, epsrel=1e-1, limit=100))

produces

(array(7.6323447), array(0.01666253))

However scipys quad method reaches an error tolerance of 1-08, which I really cant reproduce with quadpy, even when setting the limit to a really high value like 1000000

Why is that?

Update: It seems the behaviour I experienced is a bug, now reported in: https://github.com/nschloe/quadpy/issues/255

However in general this answer answers my initial question.

charelf
  • 3,103
  • 4
  • 29
  • 51