I have a very complex integral to calculate:
from __future__ import division
from scipy.integrate import quad, nquad
import numpy as np
alpha = np.array([0.298073, 1.242567, 5.782948, 38.474970])
trial = np.array([0.08704173, 0.52509737, 0.51920929, 0.31233737])
class EigenvalueProblem:
def __init__(self, a, t):
self.alpha = a
self.trial = t
# Hamiltonian, interaction part
def hartree_integrand(self, coeff):
def hartree_potential(rr2):
return np.array([coeff[ii] * coeff[jj] *
np.exp(-(self.alpha[ii] +
self.alpha[jj]) * rr2 ** 2)
for ii in range(0, 4) for jj in range(0, 4)]).sum()
def length(theta, rr1, rr2):
return 1 / np.sqrt(rr1 ** 2 + rr2 ** 2 -
2 * rr1 * rr2 * np.cos(theta))
def tmp(theta, rr1, rr2):
return 8 * np.pi ** 2 * rr1 ** 2 * rr2 ** 2 * \
np.sin(theta) * hartree_potential(rr2) * \
length(theta, rr1, rr2)
def integrand(ii, jj, theta, rr1, rr2):
return np.exp(-(self.alpha[ii] + self.alpha[jj]) * rr1 ** 2) * tmp(theta, rr1, rr2)
return [
nquad(lambda theta, rr1, rr2: integrand(i, j, theta, rr1, rr2),
[[0, np.pi], [0, np.inf], [0, np.inf]]) for i in range(0, 4) for j in range(0, 4)]
hat = EigenvalueProblem(alpha, trial)
print hat.hartree_integrand(trial)
mathematically what I want to calculate is like this (which is the integrand function), with paremeters here. However, it takes more than several hours to compute this integral. I wonder is there any method to speed up it? Thank you very much!
You ought first to extend limits in integration over r1 and r2 to be from -Infinity to +Infinity - extend limits, multiply by 1/2*1/2, etc
Second, switch to use Gauss-Hermite quadrature, which is exactly suited to integrate function with e-x2 kernels.
Appropriate code is in NumPy, see references therein
