Search code examples
pythonscipyintegral

scipy.integrate.trapz and discontinuous functions

The function scipy.integrate.trapz uses Newton-Cotes formula of order 1 as it said in the scipy documentation. However, in the derivation of this formula it is usually assumed that

  • the integrand is a continuous function and
  • the points, in which the value of the integrand is known, are distinct.

However, I tried to approximate the integral of the function f:[0,2] --> [0,2], defined by f(x) = 0 if x < 1 else 2 by calling

scipy.integrate.trapz([0, 0, 2, 2], [0, 1, 1, 2])

and obtained the right result (2.0). In the upper call,

  • an integrand is NOT a continuous function, and
  • the points in which the value of the integrand is known, are NOT distinct.

Can this "hack" be safely used in the way presented in the example?

(For each point of discontinuity x, insert x twice in the list of points, and insert the left and right limit of the integrand into the corresponding places in the list of values.)

Solution

  • The effect of repeating an x-value in trapz is the same as splitting the interval of integration at that value, and applying trapz separately to each part: 

    from scipy.integrate import trapz
trapz([0, 0, 2, 2], [0, 1, 1, 2])
trapz([0, 0], [0, 1]) + trapz([2, 2], [1, 2])

    Both results are the same. And indeed, the best way to integrate a piecewise function like yours is to split the interval of integration at the discontinuity, because then you are integrating two continuous functions. Your approach is correct.

    Theory stuff

    In some sense, your "hack" isn't necessary. The trapezoidal rule converges to the correct value of the integral for any Riemann-integrable function, which covers all piecewise continuous functions and more. That a textbook presents its proof for continuous functions only is the choice of the textbook's author, wishing to have a simpler proof. A more general theorem could be proved instead.

    That is, for any function arising in practice, the output of the trapezoidal rule will be close to its integral provided that the x-values form a sufficiently fine partition of the interval of integration. For example, uniformly spaced points will always work if you take sufficiently many.

    However, in practice one is concerned about the speed of convergence, which is noticeably worse when the function has a discontinuity. From this practical point of view, including each point of discontinuity twice, with y-values being the left and right limits, is a significant boost for the accuracy of trapezoidal rule.