How to Numerically Solve a Complex Multiple Integral in Physics?

[Korbid]

For a bidimensional system of N particles, the hamiltonian of pair-interaction is:
H(\vec{q}_1,\vec{q}_2;\vec{p}_1,\vec{p}_2)=K(\vec{p}_1,\vec{p}_2)+U(\vec{q}_1,\vec{q}_2;\vec{p}_1,\vec{p}_2)
where K is the kinetic (translational) energy and U is the potential energy
i want to solve this multiple integral:
\int\int\int\int{e^{-\frac{H(\vec{q}_1,\vec{q}_2;\vec{p}_1,\vec{p}_2)}{{k_BT}}}}d\vec{q}_1d\vec{q}_2d\vec{p}_1d\vec{p}_2
But the pair-potential depends on positions, and momentums as well:
U=\frac{k}{\tau}e^{\tau/\tau_0}
where τ₀ and κ are parameters and τ=τ(\vec{q}_{12};\vec{p}_{12})
is a function that depends on relative positions and relative momentums.
how could i solve this horrible integral? i don't need an analytical solution, a numerical solution with any software like SAGE or Mathematica is fine.

 

[jedishrfu]

Here's some Mathematica info on the topic for a start:

https://reference.wolfram.com/language/tutorial/NumericalIntegration.html

https://reference.wolfram.com/language/tutorial/NIntegrateOverview.html

and here's the equivalent info for Matlab:

http://www.mathworks.com/help/matlab/ref/integral.html

and here's some info using python and scipy:

http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html