Help integrating a vector function

[piperfw]

Hi, I'm looking to integrate a vector-valued function of a single variable. Following the example

def integrand(x):
    return [numpy.sin(x), numpy.cos(x)]

res, err = quadpy.quad(integrand, 0, 1)

I tried

def integrand(x):
    return [numpy.sin(x), 1.0]

res, err = quadpy.quad(integrand, 0, 1)

Expecting the second element of res to be 1.0, but instead I get the ValueError
ValueError: Expected the function return value to be of shape (..., 21), but got shape (2,) instead.

Where am I going wrong? Apologies for such a basic question. I noted that changing 1.0 to 1.0+x*0.0 allowed the second integral to be performed.

[nschloe]

import numpy as np
import quadpy


def integrand(x):
    return [np.sin(x), np.ones_like(x)]

res, err = quadpy.quad(integrand, 0, 1)

[nschloe]

With imports, please. I just want to copy and paste and note waste time guessing the modules.

[nschloe]

So, t is of shape (21,), i.e., you're asked to return your function evaluated at 21 points. You can do that via a loop to start with. Later, if you're about speed, you want to vectorize it. This

    out = np.array([comm(f(T)) for T in t])

has shape (21, 4, 4), but you need shape (..., 21). One solution is to just move the first axis to the end:

def g(t):
    out = np.array([comm(f(T)) for T in t])
    return np.moveaxis(out, 0, -1)

After vectorization, all these shenanigans will disappear and your code will be readable and fast.
Anyway, the result is

[[ 0. -1.  1.  0.]
 [-1.  2.  0.  1.]
 [ 1.  0. -2. -1.]
 [ 0.  1. -1.  0.]]

[piperfw]

@nschloe that is really helpful, thank you for working on this problem with me. May I check - the number '21' isn't magic but just some minimum number of points the integration scheme requires in order to ensure small error? So in more difficult cases one would expect it to be larger.

[nschloe]

21 Is the default number of points in the integration scheme; it's the 21-point Gauss–Kronrod quadrature. Because this integral here is so trivial, it doesn't need to subdivide the interval. If the function is more complex, expect many more points.

[piperfw]

That makes sense. Let me know if there is any simple docs/wiki I might be able to contribute to.