Numerical calculation of triple integral using Python

Question

I am in deep need to calculate this integral. I have been trying to do so for a couple of months, using the numpy package in Python, and in particular, the integrate.tplquad function.

from __future__ import division
from math import *
import numpy as np
import scipy.special as special
import scipy.integrate as integrate

a=1.e-19
b=1.e-09
zo=1.e7
H=1.e15

v=1.e18

def integrand(v,z,x,u):
    value=x**(-0.5)*special.kv(5/3,u)*(a*v*z/x-1/2.)*exp(-b*sqrt(z*v/x))
    return value

i=integrate.tplquad(lambda u,x,z: integrand(v,z,x,u),1.e7,1.e15,lambda z:0.,lambda z:np.inf, lambda x,z : x, lambda x,z : np.inf)

print i

In the code above, I have tried a value of v=10^18, in order to normalize the argument of the exponent, and not get too small or too large coefficients.

However, no matter what value of v I plug in, I always get

out: (0.0, 0.0)

I do not know how to surpass this problem.

I have also tried expanding the exponential function into a power series, but I get the same result.

Now, I know for a fact that the integral must have a finite, positive value for all v. In fact, I'd be delighted if I could calculate it for any v.

If anyone has encountered a similar problem, I'd be delighted if they could share their wisdom. Any help is welcome. Thanks.

Answer 1

First thing is that the u and z integrals can be solved exactly. The result is a rather convoluted function involving exponentials, the gamma function, and generalized hyper geometric series. The advantage is that it's only on one variable, and so it can be easily examined graphically. Here are some of the curves, for different values of \nu:

And here is the expression:

It's convenient to integrate this function, as it's much faster and exact to do so. But, and here's the second point, this suffers from numerical issues due to machine precision as x -> \inf. Here are a couple of plots clearly showing the issues:

When plotting instead with arbitrary working precision, the problem disappears:

So, the numerical issues have to be dealt with too, by using an arbitrary precision library like mpmath under Python, and/or by ignoring/discarding the upper leg of the integration interval, i.e. in this case by example, integrating between 0 and 19 / 20, instead of 0 and \inf.

Below is a Python program which, using mpmath, integrates the expression above(an equivalent one) between x=0 and x=20

#!/usr/bin/env python3
#encoding: utf

from mpmath import mp, mpf, sqrt, besselk, exp, quad, pi, hyper, gamma

maxprecision = 64 # significant digits
maxdegree    = 3  # maximum degree of the quadrature rule

mp.dps = maxprecision

# z0 = mpf(1.e7)
# H  = mpf(1.e15)
a  = mpf(1.e-19)
b  = mpf(1.e-9)

sqrt3 = sqrt(3.)
sqrt10 = sqrt(10.)
inf = mpf('inf')

epsilon=10.**-maxprecision

def integrand(z, x, u):
  value = 1./sqrt(x) * besselk(5./3, u) * (a*z*nu/x - 1./2) * exp(-b * sqrt(z*nu/x))
  return value

def integrand3(x):
  value = 1. / (960. * b**4 * x**(19./6) * (nu / x)**(3./2)) * exp(-10000000. * sqrt10 * b * sqrt(nu/x)) * (-b**2 * x * (1000. * sqrt10 * b * (-10000. + exp(9999000. * sqrt10 * b * sqrt(nu/x))) * nu + (-1. + exp(9999000. * sqrt10 * b * sqrt(nu/x))) * x * sqrt(nu/x)) + 4. * a * (3000. * sqrt10 * b * (-10000. + exp(9999000. * sqrt10 * b * sqrt(nu/x))) * x * nu + 5000000000. * sqrt10 * b**3 * (-1000000000000. + exp(9999000. * sqrt10 * b * sqrt(nu/x))) * nu**2 + 3. * (-1. + exp(9999000. * sqrt10 * b * sqrt(nu/x))) * x**2 * sqrt(nu/x) + 15000000. * b**2 * (-100000000. + exp(9999000. * sqrt10 * b * sqrt(nu/x))) * x * nu * sqrt(nu/x))) * (-320. * sqrt3 * pi * x**(2./3) + 960. * 2.**(2./3) * gamma(2./3) * hyper([-1./3], [-2./3, 2./3], x**2 / 4.) + 27. * 2.**(1./3) * x**(10./3) * gamma(-2./3) * hyper([4./3], [7./3, 8./3], x**2 / 4.))
  return value

for e in range(0, 19):
  nu = mpf(10**e)
#   I1 = quad(lambda x: quad(lambda u, z: integrand(z, x, u), [x, inf], [z0, H], method='tanh-sinh', maxdegree=maxdegree), [0., inf], method='tanh-sinh', maxdegree=maxdegree)
#   print("ν = 10^%d: NI(x, u, z) = %f" % (e, I1))
  I3 = quad(lambda x: integrand3(x), [0., 20.], method='tanh-sinh', maxdegree=maxdegree)
  print("ν = 10^%d: NI(x)       = %f" % (e, I3))
#   print("ν = 10^%d: error       = %.2f%% " % (e, (I3-I1)/(I1+epsilon)*100.))

The results are:

ν = 10^0: NI(x)       = -12118569382494810.000000
ν = 10^1: NI(x)       = -6061688705992958.000000
ν = 10^2: NI(x)       = -2359248732202052.500000
ν = 10^3: NI(x)       = -535994574128436.812500
ν = 10^4: NI(x)       = -26417279314541.281250
ν = 10^5: NI(x)       = 3636613281577.158203
ν = 10^6: NI(x)       = 416805025513.477356
ν = 10^7: NI(x)       = 41860949922.522430
ν = 10^8: NI(x)       = 4285965504.873075
ν = 10^9: NI(x)       = 477094892.498829
ν = 10^10: NI(x)       = 65240304.226613
ν = 10^11: NI(x)       = 9524738.222360
ν = 10^12: NI(x)       = 680659.198974
ν = 10^13: NI(x)       = 5287.165984
ν = 10^14: NI(x)       = 0.224778
ν = 10^15: NI(x)       = 0.000000
ν = 10^16: NI(x)       = -0.000000
ν = 10^17: NI(x)       = -0.000000
ν = 10^18: NI(x)       = -0.000000