2016-09-03 2 views
1

Comment évaluer xe^x/(e^x-1) avec une stabilité numérique autour de zéro et quand x est très positif ou négatif? J'ai accès à toutes les fonctions mathématiques habituelles dans numpy et scipy.Comment évaluer xe^x/(e^x-1) avec la stabilité numérique en Python?

+0

Juste trouvé scipy.special.exprel, qui fonctionnerait. J'ai besoin que mon code fonctionne sur theano ou tensorflow aussi, donc une solution plus basique serait meilleure. –

+0

Frustrant, 'np.expm1' semble donner' nan' au lieu de 'inf' pour les grandes entrées. Cela rend plus difficile d'éviter des choses comme «si» ou «où». – user2357112

+0

Cela semble être un [bug dépendant de la plateforme] (https://github.com/numpy/numpy/issues/6818), toujours ouvert. – user2357112

Répondre

2
def f(x): 
    if abs(x) > 0.1: return x*exp(x)/(exp(x)-1) 
    else: return 1/(1.-x/2.+x**2/6.-x**3/24.) 

L'expansion dans la dernière ligne peut être prolongée de la façon évidente si une plus grande précision est nécessaire, et peut être plus rapidement en re-phrasé. En l'état, il se trompe autant que 1e-6.

+2

Il y a une valeur pas si grande 'z' telle que pour' x> z' la valeur réelle de l'expression 'x * exp (x)/(exp (x) -1)' est très proche de 'x', alors que le calculer directement entraînera un débordement et nan. – Leon

+0

Bon point, @Leon. Pour gérer cette possibilité, les cas de 'x> 0.1' et' x <-0.1' devraient être traités séparément. Si 'x> 0.1', retourner' x/(1-exp (-x)) '; si 'x <-0.1', retourner' x * exp (x)/(exp (x) -1) '; et sinon utiliser l'approximation polynomiale. –