Calcul numérique de l'intégrale triple en utilisant Python

Question

Je suis dans le besoin profond de calculer this intégrale. J'ai essayé de le faire pendant quelques mois, en utilisant le paquet numpy en Python, et en particulier, la fonction 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 

Dans le code ci-dessus, j'ai essayé une valeur de v = 10^18, afin de normaliser l'argument de l'exposant, et ne pas obtenir les coefficients trop petits ou trop grands.

Cependant, quelle que soit la valeur de v je branche, je reçois toujours

out: (0.0, 0.0) 

Je ne sais pas comment dépasser ce problème.

J'ai également essayé d'étendre la fonction exponentielle dans une série de puissance, mais j'ai le même résultat.

Maintenant, je sais pertinemment que l'intégrale doit avoir une valeur finie, positive pour tous v. En fait, je serais heureux si je pouvais calculer pour tout v.

Si quelqu'un a rencontré un problème similaire, je serais ravi s'ils pouvaient partager leur sagesse. Toute aide est la bienvenue. Merci.

Answer 1

La première chose à faire est que les intégrales u et z peuvent être résolues exactement. Le résultat est une fonction plutôt alambiquée impliquant des exponentielles, la fonction gamma et des séries hyper géométriques généralisées. L'avantage est que c'est seulement sur une variable, et donc il peut être facilement examiné graphiquement. Voici quelques-unes des courbes, pour différentes valeurs de \ nu:

Et voici l'expression:

Il est pratique d'intégrer cette fonction, car il est beaucoup plus rapide et exact de le faire. Mais, et voici le deuxième point, cela souffre de problèmes numériques dus à la précision de la machine comme x -> \ inf. Voici quelques parcelles montrant clairement les enjeux:

Lors du traçage au lieu avec une précision de travail arbitraire, le problème disparait:

Ainsi, le numérique Les problèmes doivent également être résolus en utilisant une bibliothèque de précision arbitraire comme mpmath sous Python, et/ou en ignorant/en supprimant t la partie supérieure de l'intervalle d'intégration, c'est-à-dire dans ce cas par exemple, intégrant entre 0 et 19/20, au lieu de 0 et \ inf.

est Ci-dessous un programme Python qui, en utilisant mpmath, intègre l'expression ci-dessus (un poste équivalent) entre x = 0 et 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.)) 

Les résultats sont les suivants:

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