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の値を試してみました。しかし

、関係なく、私はいつも、私はこの問題を凌駕する方法がわからない

out: (0.0, 0.0) 

を私はプラグインVのどのような値を取得していません。

私は指数関数をべき級数に拡張しようとしましたが、同じ結果が得られます。

さて、私は積分はすべてVのための有限正の値を持たなければならないという事実を知っている。私はVのためにそれを計算することができれば実際に、私は喜んだろう。

誰もが遭遇した場合同様の問題で、彼らの知恵を分かち合うことができれば嬉しいです。どんな助けも歓迎です。ありがとう。

Answer 1

最初に、uとzの積分を正確に解くことができます。その結果は、指数関数、ガンマ関数、および一般化された超幾何系列を含む複雑な関数である。利点は、それが1つの変数にしかないことであり、したがって、簡単にグラフィカルに検査することができるということです。ここ

と表現されています：

それはそれがはるかだとして、この機能を統合すると便利です。ここの\ NUの異なる値の曲線の一部であり、迅速かつ正確に行うことができます。しかし、ここでは2番目の点がありますが、これはx - > \ infのような機械精度による数値的問題に苦しんでいます。

任意の加工精度を代わりにプロットし、問題が消える：

ので、数値ここで明確な問題を示すプロットのカップルですPythonの下でmpmathのような任意の精度ライブラリを使用すること、および/またはtを無視/破棄することによっても問題を処理する必要があります積分区間の上肢、すなわち、この例では、0と\ infの代わりに、0と19/20との間に積分する。以下

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.)) 

間の結果であるmpmathを使用して、（同等のもの）上記式を積分し、Pythonプログラムである：

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