numpyの指数関数のオーバーフロー

Question

実現したいこと

実験データをまとめるために、
ある程度大きな数でも、指数関数の計算を行いたい。

前提

実験データをまとめるために、
指数関数を用いた計算を行なうコードを作成していたところ、
24行目のnp.exp()がオーバーフローしてしまった。
(コード中のデータは、オンライン上で一般に公開されているデータを使用しています。)

• 1本目はオーバーフロー以下のため、計算できた。
• 2本目〜5本目はオーバーフローが発生し、計算できなかった。

2本目〜5本目の計算に耐えうる計算の方法が知りたいです。

発生している問題・エラーメッセージ

RuntimeWarning: overflow encountered in exp
  denom = np.exp(const_b * const_c * (x ** (-1) - 1)) - 1

import numpy as np
import matplotlib.pyplot as plt

#x：濃度
#xの範囲を設定
x = np.arange(0.01, 0.07, 0.001)

#定数
#Aτ0＝const_a
#B=const_b

#c:密度
#c=const_c
const_c = 2.6


#y:τy
#分子:numer
#分母:denom

#yを計算する関数を定義
def stress_y(x, const_a, const_b, const_c):
    numer = (const_c * (x ** (-1)) + (1 - const_c)) ** (-2)
    denom = np.exp(const_b * const_c * (x ** (-1) - 1)) - 1
    y = const_a * (numer / denom)
    return y

#濃度に応じた値を計算
#1本目(0.011M)
const_a1 = 87.83
const_b1 = 7.03

y1 = stress_y(x, const_a1, const_b1, const_c)

#2本目(0.022M)
const_a2 = 741.71
const_b2 = 26.96

y2 = stress_y(x, const_a2, const_b2, const_c)

#3本目(0.033M)
const_a3 = 3617.37
const_b3 = 55.46

y3 = stress_y(x, const_a3, const_b3, const_c)

#4本目(0.044M)
const_a4 = 16338.03
const_b4 = 88.04

y4 = stress_y(x, const_a4, const_b4, const_c)

#5本目(0.055M)
const_a5 = 302752.7
const_b5 = 173.25

y5 = stress_y(x, const_a5, const_b5, const_c)

ソースコード

試したこと

対数関数の使用は思いついたのですが、具体的な実装方法がわかりません。

補足情報（FW/ツールのバージョンなど）

spyderを使用しています。

Comments

[melian]

dtype に np.float128 を指定すると 2本目はオーバーフローしないのですが、3〜5本目はオーバーフローしてしまいますね…

x = np.arange(0.01, 0.07, 0.001, dtype=np.float128)

[8524ba23]

- その式で何を計算しているのか
- 各変数の意味と単位
を質問本文に追記すると回答得られやすくなると思います。
また、1本目の計算で得られた結果は正しかったということでしょうか？

[bsdfan]

計算結果はプロットしたいのでしょうか？
最終的にどういうアウトプットが欲しいかによって手段を選ぶ必要があります。

Answer 1

mpmath - Python library for arbitrary-precision floating-point arithmetic を利用する場合。

import numpy as np
import mpmath

mpmath.mp.dps = 50

#x：濃度
#xの範囲を設定
x = np.arange(0.01, 0.07, 0.001)

#定数
const_c = 2.6

#yを計算する関数を定義
power_array = np.frompyfunc(mpmath.power, 2, 1)
exp_array = np.frompyfunc(mpmath.exp, 1, 1)
def stress_y(x, const_a, const_b, const_c):
    numer = power_array(const_c * power_array(x, -1) + (1 - const_c), -2)
    denom = exp_array(const_b * const_c * (power_array(x, -1) - 1)) - 1
    y = const_a * (numer / denom)
    return y

#濃度に応じた値を計算
#1本目(0.011M)

        :

# [mpf('1.7932398509323452135403010297111599472201279379812348e-789')
#  mpf('3.1689260903784489600004401703364548794600593787163488e-717')
#  mpf('5.1719326727630832701870953164567431379194547158702907e-657')
#                      :
#