python - 使用Python进行三重积分的数值计算

Question

我很需要计算this不可缺少的。几个月来，我一直在尝试这样做，使用 Python 中的 numpy 包，特别是 integrate.tplquad 函数。

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

在上面的代码中，我尝试了一个 v=10^18 的值，以便对指数的参数进行归一化，并且不会得到太小或太大的系数。

但是，无论我插入什么值的v，我总是得到

out: (0.0, 0.0)

我不知道如何超越这个问题。

我也试过将指数函数展开为幂级数，但我得到了相同的结果。

现在，我知道一个事实，即所有 v 的积分必须有一个有限的正值。事实上，如果我能计算任何 v 的积分，我会很高兴。

如果有人遇到过类似的问题，很高兴能分享他们的智慧。欢迎任何帮助。谢谢。

Answer 1

首先是可以精确求解 u 和 z 积分。结果是一个相当复杂的函数，涉及指数、 Gamma 函数和广义超几何级数。优点是它只针对一个变量，因此可以很容易地以图形方式对其进行检查。下面是一些曲线，对于不同的\nu 值:

下面是表达式:

集成这个函数很方便，因为这样做更快更准确。但是，这是第二点，由于 x ->\inf 的机器精度，这会遇到数值问题。这里有几个图表清楚地显示了问题:

当以任意工作精度绘图时，问题消失了:

因此，数值问题也必须通过在 Python 下使用任意精度库(如 mpmath)和/或通过忽略/丢弃积分区间的上段来处理，即在这个例子，在 0 和 19/20 之间积分，而不是 0 和\inf。

下面是一个 Python 程序，它使用 mpmath 在 x=0 和 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.))

结果是:

ν = 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