使用Python進行三元積分的數值計算

我深深需要計算this積分。我一直試圖這樣做幾個月，使用Python中的numpy包，特別是integrate.tplquad函數。使用Python進行三元積分的數值計算

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有限，正值。事實上，我會很高興，如果我能計算出它的任何訴

如果有人遇到一個類似的問題，如果他們能分享他們的智慧，我會很高興。歡迎任何幫助。謝謝。

來源

2016-11-22 Valentina

你應該在腳本的頂部有'from __future__ import division'。否則'special.kv'中的'5/3'是'1'，而不是'1.666667'。（儘管如此，這並不能解決整個問題）。 – Elliot

哦對，謝謝你注意:) – Valentina

你的積分在x = 0處有一個奇點，這是積分的限制之一。 –

首先，u和z積分可以完全解決。結果是一個相當複雜的函數，涉及指數，伽馬函數和廣義超幾何級數。優點是它只在一個變量上，因此可以很容易地用圖形方式進行檢查。以下是一些曲線，爲\ NU的不同值：

，這裏是表達：

它的方便集成這個功能，因爲它是多更快，更精確地這樣做。但是，這裏是第二點，由於機器精度而產生數值問題，如x - > \ inf。這裏有幾個地塊清晰顯示問題：

當任意加工精度，而不是陰謀，問題就會消失：

因此，數值問題也必須處理，通過在Python下使用像mpmath這樣的任意精度庫，和/或通過忽略/丟棄t積分區間的上限，即在這種情況下通過示例，在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

來源

2016-12-01 16:49:42

使用Python進行三元積分的數值計算

回答

相關問題