[Python] Why does my time-dependent integration using Cubature in Python produce non-constant...

I'm trying to numerically evaluate a multidimensional integral using the cubature package in Python. My goal is to compute a rate function that should be constant in time, but the output varies with time.

Here’s a summary of what I’m doing:

The integral has three levels:

Innermost: angular integration over θ and φ

Middle: radial integration over momentum q

Outermost: energy integration over ω

The integrand includes a velocity term vlab(t) that changes with time, but the physical quantity (rate) is expected to be constant due to isotropy.

The function g(...), which depends on vlab(t), includes a sharp cutoff (vminus > vesc) that might be introducing numerical instability.

I'm using nested calls to cubature. Despite the function being theoretically time-independent after integration, I'm observing time-varying results. Here's the relevant part of my code:

def g(qunit, q, w, t):
vminus = w/q + q/(2*mx) + np.dot(qunit, vlab(t))
if vminus > vesc:
return 0
return (np.pi * vo**2 / (q * No)) * (np.exp(-vminus**2 / vo**2) - np.exp(-vesc**2 / vo**2))

def integrand(q2d, q, w, t):
theta = q2d[0]
phi = q2d[1]
qvec = np.array([q*np.sin(theta)*np.cos(phi), q*np.sin(theta)*np.sin(phi),
q*np.cos(theta)])
qunit = qvec/q
jacobian = q**2*np.sin(theta)
val = sigma * g(qunit, q, w, t) * jacobian * somefunction(q)
return val



What I’ve tried: Increasing maxEval and decreasing relerr and abserr.

Continue reading...