27

Motivation: I have a multidimensional integral, which for completeness I have reproduced below. It comes from the computation of the second virial coefficient when there is significant anisotropy:

enter image description here

Here W is a function of all the variables. It is a known function, one which I can define a python function for.

Programming Question: How do I get scipy to integrate this expression? I was thinking of chaining two triple quads (scipy.integrate.tplquad) together, but I'm worried about performance and accuracy. Is there a higher dimensional integrator in scipy, one that can handle an arbitrary number of nested integrals? If not, what is the best way to do this?

Improve this question
1
  • You may be better try Sympy. Commented Dec 28, 2012 at 16:03

4 Answers 4

Reset to default
36

With a higher-dimensional integral like this, monte carlo methods are often a useful technique - they converge on the answer as the inverse square root of the number of function evaluations, which is better for higher dimension then you'll generally get out of even fairly sophisticated adaptive methods (unless you know something very specific about your integrand - symmetries that can be exploited, etc.)

The mcint package performs a monte carlo integration: running with a non-trivial W that is nonetheless integrable so we know the answer we get (note that I've truncated r to be from [0,1); you'll have to do some sort of log transform or something to get that semi-unbounded domain into something tractable for most numerical integrators):

import mcint import random import math def w(r, theta, phi, alpha, beta, gamma): return(-math.log(theta * beta)) def integrand(x): r = x[0] theta = x[1] alpha = x[2] beta = x[3] gamma = x[4] phi = x[5] k = 1. T = 1. ww = w(r, theta, phi, alpha, beta, gamma) return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta) def sampler(): while True: r = random.uniform(0.,1.) theta = random.uniform(0.,2.*math.pi) alpha = random.uniform(0.,2.*math.pi) beta = random.uniform(0.,2.*math.pi) gamma = random.uniform(0.,2.*math.pi) phi = random.uniform(0.,math.pi) yield (r, theta, alpha, beta, gamma, phi) domainsize = math.pow(2*math.pi,4)*math.pi*1 expected = 16*math.pow(math.pi,5)/3. for nmc in [1000, 10000, 100000, 1000000, 10000000, 100000000]: random.seed(1) result, error = mcint.integrate(integrand, sampler(), measure=domainsize, n=nmc) diff = abs(result - expected) print "Using n = ", nmc print "Result = ", result, "estimated error = ", error print "Known result = ", expected, " error = ", diff, " = ", 100.*diff/expected, "%" print " "

Running gives

Using n = 1000 Result = 1654.19633236 estimated error = 399.360391622 Known result = 1632.10498552 error = 22.0913468345 = 1.35354937522 % Using n = 10000 Result = 1634.88583778 estimated error = 128.824988953 Known result = 1632.10498552 error = 2.78085225405 = 0.170384397984 % Using n = 100000 Result = 1646.72936 estimated error = 41.3384733174 Known result = 1632.10498552 error = 14.6243744747 = 0.8960437352 % Using n = 1000000 Result = 1640.67189792 estimated error = 13.0282663003 Known result = 1632.10498552 error = 8.56691239895 = 0.524899591322 % Using n = 10000000 Result = 1635.52135088 estimated error = 4.12131562436 Known result = 1632.10498552 error = 3.41636536248 = 0.209322647304 % Using n = 100000000 Result = 1631.5982799 estimated error = 1.30214644297 Known result = 1632.10498552 error = 0.506705620147 = 0.0310461413109 %

You could greatly speed this up by vectorizing the random number generation, etc.

Of course, you can chain the triple integrals as you suggest:

import numpy import scipy.integrate import math def w(r, theta, phi, alpha, beta, gamma): return(-math.log(theta * beta)) def integrand(phi, alpha, gamma, r, theta, beta): ww = w(r, theta, phi, alpha, beta, gamma) k = 1. T = 1. return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta) # limits of integration def zero(x, y=0): return 0. def one(x, y=0): return 1. def pi(x, y=0): return math.pi def twopi(x, y=0): return 2.*math.pi # integrate over phi [0, Pi), alpha [0, 2 Pi), gamma [0, 2 Pi) def secondIntegrals(r, theta, beta): res, err = scipy.integrate.tplquad(integrand, 0., 2.*math.pi, zero, twopi, zero, pi, args=(r, theta, beta)) return res # integrate over r [0, 1), beta [0, 2 Pi), theta [0, 2 Pi) def integral(): return scipy.integrate.tplquad(secondIntegrals, 0., 2.*math.pi, zero, twopi, zero, one) expected = 16*math.pow(math.pi,5)/3. result, err = integral() diff = abs(result - expected) print "Result = ", result, " estimated error = ", err print "Known result = ", expected, " error = ", diff, " = ", 100.*diff/expected, "%"

which is slow but gives very good results for this simple case. Which is better is going to come down to how complicated your W is and what your accuracy requirements are. Simple (fast to evaluate) W with high accuracy will push you to this sort of method; complicated (slow to evaluate) W with moderate accuracy requirements will push you towards MC techniques.

Result = 1632.10498552 estimated error = 3.59054059995e-11 Known result = 1632.10498552 error = 4.54747350886e-13 = 2.7862628625e-14 %
Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

Thanks! I'll take a look at mcint and see if it performs better than my ad-hoc MC method I've got going now.
@JohnathanDursi is it possible to get a multidimensional Gaussian quadrature in Python? Such quadrature sets are used for instance in solving the heat conduction equation. In this case one distributes the polar angles according to some quadrature rule and the azimuthal angles (directions) are uniformly distributed.
14

Jonathan Dursi has made a very good answer. I will just add on to his answer.

Now scipy.integrate has a function named nquad that one can perform a multi-dimensional integral without hassle. See this link for more information. Below we compute the integral using nquad with Jonathan's example:

from scipy import integrate import math import numpy as np def w(r, theta, phi, alpha, beta, gamma): return(-math.log(theta * beta)) def integrand(r, theta, phi, alpha, beta, gamma): ww = w(r, theta, phi, alpha, beta, gamma) k = 1. T = 1. return (math.exp(-ww/(k*T)) - 1.)*r*r*math.sin(beta)*math.sin(theta) result, error = integrate.nquad(integrand, [[0, 1], # r [0, 2 * math.pi], # theta [0, math.pi], # phi [0, 2 * math.pi], # alpha [0, 2 * math.pi], # beta [0, 2 * math.pi]]) # gamma expected = 16*math.pow(math.pi,5)/3 diff = abs(result - expected)

The result is more accurate than the chained tplquad:

>>> print(diff) 0.0
Improve this answer

Comments

10

I'll just make a couple of general comments about how to accurately do this sort of integral, but this advice is not specific to scipy (too long for a comment, even though it is not an answer).

I don't know your use case, i.e. whether you are satisfied with a `good' answer with a few digits of accuracy which could be obtained straightforwardly using Monte Carlo as outlined in Jonathan Dursi's answer, or whether you really want to push the numerical accuracy as far as possible.

I've performed analytic, Monte Carlo and quadrature calculations of virial coefficients myself. If you want to do the integrals accurately, then there are a few things you should do:

  1. Attempt to perform as many of the integrals exactly as possible; it may well be that integration in some of your coordinates is quite simple.

  2. Consider transforming your variables of integration so that the integrand is as smooth as possible. (This helps for both Monte Carlo and quadrature).

  3. For Monte Carlo, use importance sampling for best convergence.

  4. For quadrature, with 7 integrals it may just be possible to get really fast convergence using tanh-sinh quadrature. If you can get it down to 5 integrals then you should be able to get 10's of digits of precision for your integral. I highly recommend mathtool / ARPREC for this purpose, available from David Bailey's homepage: http://www.davidhbailey.com/

Improve this answer

1 Comment

Thanks for the input. Do you mind on elaborating on #2? A priori how do I know what a good transform would be? Since you've done these type of calculations before, any additional input would be appreciated.
1

First to say that I am not that good in math so please be kind. Anyway, here is my try:
Note that in your question there are 6 variables but 7 integrals!?
In Python using Sympy: 

>>> r,theta,phi,alpha,beta,gamma,W,k,T = symbols('r theta phi alpha beta gamma W k T') >>> W = r+theta+phi+alpha+beta+gamma >>> Integral((exp(-W/(k*T))-1)*r**2*sin(beta)*sin(theta),(r,(0,2*pi)),(theta,(0,pi)),(phi,(0,2*pi)),(alpha,(0,2*pi)),(beta,(0,pi)),(gamma,(0,pi))) >>> integrate((exp(-W)-1)*r**2*sin(beta)*sin(theta),(r,(0,2*pi)),(theta,(0,pi)),(phi,(0,2*pi)),(alpha,(0,2*pi)),(beta,(0,pi)),(gamma,(0,pi)))

and here is the result: [LateX code]

\begin{equation*}- \frac{128}{3} \pi^{6} - \frac{\pi^{2}}{e^{2 \pi}} - \frac{\pi}{e^{2 \pi}} - \frac{2}{e^{2 \pi}} - \frac{\pi^{2}}{e^{3 \pi}} - \frac{\pi}{e^{3 \pi}} - \frac{2}{e^{3 \pi}} - 3 \frac{\pi^{2}}{e^{6 \pi}} - 3 \frac{\pi}{e^{6 \pi}} - \frac{2}{e^{6 \pi}} - 3 \frac{\pi^{2}}{e^{7 \pi}} - 3 \frac{\pi}{e^{7 \pi}} - \frac{2}{e^{7 \pi}} + \frac{1}{2 e^{9 \pi}} + \frac{\pi}{e^{9 \pi}} + \frac{\pi^{2}}{e^{9 \pi}} + \frac{1}{2 e^{8 \pi}} + \frac{\pi}{e^{8 \pi}} + \frac{\pi^{2}}{e^{8 \pi}} + \frac{3}{e^{5 \pi}} + 3 \frac{\pi}{e^{5 \pi}} + 3 \frac{\pi^{2}}{e^{5 \pi}} + \frac{3}{e^{4 \pi}} + 3 \frac{\pi}{e^{4 \pi}} + 3 \frac{\pi^{2}}{e^{4 \pi}} + \frac{1}{2 e^{\pi}} + \frac{1}{2}\end{equation*}

You may play a bit more for your question ;)

Improve this answer

1 Comment

That still looks like it is doing symbolic calculation, i.e. your W is a linear function of the input variables, hence the exact result. For me W is non-linear and not expressible as a math function, but as the result of another calculation (thus defined as a python function). You are correct that I should only have 6 integrals, I must have gotten carried away TeXing it up.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.