python求三重积分_三重积分的Python数值计算

首先，u积分和z积分可以精确求解。结果是一个相当复杂的函数，包括指数、伽马函数和广义超几何级数。它的优点是它只在一个变量上，因此可以很容易地用图形来检查。以下是一些不同值的\nu曲线：

下面是一个表达式：

集成这个函数很方便，因为这样做更快、更精确。但是，这是第二点，由于机器精度如x-&gt；\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.))

结果如下：

^{pr2}$