Pythonでの積分について、integrate.quadを用いて-π/2からπ/2の範囲での積分をしたい

Question

実現したいこと

Pythonでの積分について、integrate.quadを用いて-π/2からπ/2の範囲での積分をしたい

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

●RuntimeWarning: invalid value encountered in double_scalars
y_in = lambda phi: np.sin(phi)/(phi*(delta_in+Do*(1-np.cos(phi))))

●IntegrationWarning: The occurrence of roundoff error is detected, which prevents
the requested tolerance from being achieved. The error may be
underestimated.
integ_in = integrate.quad(y_in, (-1/2)*np.pi, (1/2)*np.pi)

●RuntimeWarning: invalid value encountered in double_scalars
y_out = lambda phi: np.sin(phi)/(phi*(delta_out+Do*(1-np.cos(phi))))

●IntegrationWarning: The occurrence of roundoff error is detected, which prevents
the requested tolerance from being achieved. The error may be
underestimated.
integ_out = integrate.quad(y_out, (-1/2)*np.pi, (1/2)*np.pi)

該当のソースコード

import numpy as np
from scipy import integrate

# 変数を定義
delta_in = 0.15E-3
delta_out = 1.00E-3
delta_tc = 0.6E-3
Do = 1.03E-3
Di = 1.00E-3
eps_0 = 8.854187E-12
eps_r = 4.0
h_core = 12.3E-3
dia_in_core = 17.6E-3
dia_out_core = 32.8E-3

r_in = (dia_in_core/2) - delta_tc - Do/2
r_out = (dia_out_core/2) + delta_tc + Do/2


# 積分１
y_in = lambda phi: np.sin(phi)/(phi*(delta_in+Do*(1-np.cos(phi))))
y_out = lambda phi: np.sin(phi)/(phi*(delta_out+Do*(1-np.cos(phi))))
integ_in = integrate.quad(y_in, (-1/2)*np.pi, (1/2)*np.pi)
integ_out = integrate.quad(y_out, (-1/2)*np.pi, (1/2)*np.pi)

L_T_in = h_core + delta_tc + Do
L_T_out = L_T_in
L_T = h_core*2+ (dia_out_core-dia_in_core) + Do



C_tt_air_in = (Do/2.0)*eps_0*L_T_in * integ_in[0]
C_tt_air_out = (Do/2.0)*eps_0*L_T_out * integ_out[0]



# 積分２
f = lambda phi, r: np.sin(phi)/(phi*( ((delta_out-delta_in)/(r_out-r_in))*(r-r_in + delta_in + Do*(1-np.cos(phi)))))
val, err = integrate.dblquad(f, r_in, r_out, lambda phi: (-1/2)*np.pi, lambda phi: (1/2)*np.pi)

C_tt_air_l = (Do/2.0)*eps_0 * val
print("C_tt_air_l", C_tt_air_l)

C_tt_ins = (np.pi*eps_0*eps_r*L_T)/(2*np.log(Do/Di))
print("C_tt_ins", C_tt_ins)

C_tt_in = (C_tt_ins*C_tt_air_in)/(C_tt_ins+C_tt_air_in)
C_tt_out = (C_tt_ins*C_tt_air_out)/(C_tt_ins + C_tt_air_out)
C_tt_l = (C_tt_ins*C_tt_air_l)/(C_tt_ins + C_tt_air_l)

C_tt = C_tt_in + C_tt_out + 2*C_tt_l

試したこと

積分範囲を、0からnp.pi、(-1/2)*np.piから0などにすると積分できましたが、
(-1/2)*np.piから(1/2)*np.piに関しては上記のようなエラーが出ます。

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

コードを全文掲載しました。

コード中の式は、工学系論文の式を引用しています。
引用した論文では、最終的にC_ttの値は、上記の変数の時におよそ1.124E-12の結果となるようです。

Comments

[jeanbiego]

<code>機能を使って、コードを囲ってください。読みやすくなります。
また、気になる部分の抜粋ではなくて、コピペだけでテストできるように、なるべく全体を記載したほうが有用な回答がつきやすいと思います。

[mag0123]

jeanbiego様
修正依頼ありがとうございました。コードの掲載について修正しました。

Answer 1

pythonの計算精度の問題で、途中で計算しきれなくなっているようです。（参考：15. 浮動小数点演算、その問題と制限）sin(φ)/φをsinc(φ/π)と置き換え（参考：numpy.sinc）て計算を省力化したところ、最後まで回るようになりました。（最初のコードでは、integ_inなどが計算できずにNaNになっていました）
しかし結果は、1.6141283310668653e-12ということで、まだ目標とは50%くらいずれていますね。途中の誤差の積み重ねでしょうか。精度を上げるには、途中の式を手でほぐして簡単な計算に置き換えていくのが王道かと思います。

import numpy as np
from scipy import integrate

# 変数を定義
delta_in = 0.15E-3
delta_out = 1.00E-3
delta_tc = 0.6E-3
Do = 1.03E-3
Di = 1.00E-3
eps_0 = 8.854187E-12
eps_r = 4.0
h_core = 12.3E-3
dia_in_core = 17.6E-3
dia_out_core = 32.8E-3

r_in = (dia_in_core/2) - delta_tc - Do/2
r_out = (dia_out_core/2) + delta_tc + Do/2
# r_in: 0.007685000000000001
# r_out: 0.017515000000000003

# 積分１
y_in = lambda phi: np.sinc(phi/np.pi)/(delta_in+Do*(1-np.cos(phi)))
y_out = lambda phi: np.sinc(phi/np.pi)/(delta_out+Do*(1-np.cos(phi)))

integ_in = integrate.quad(y_in, (-1/2)*np.pi, (1/2)*np.pi)
integ_out = integrate.quad(y_out, (-1/2)*np.pi, (1/2)*np.pi)
# integ_in: (8546.321590461577, 9.466839635601456e-05)
# integ_out: (2154.5245517896133, 2.3920027417134726e-11)


L_T_in = h_core + delta_tc + Do
L_T_out = L_T_in
L_T = h_core*2+ (dia_out_core-dia_in_core) + Do
# L_T_in: 0.01393
# L_T_out: 0.01393

C_tt_air_in = (Do/2.0)*eps_0*L_T_in * integ_in[0]
C_tt_air_out = (Do/2.0)*eps_0*L_T_out * integ_out[0]
# C_tt_air_in: 5.428580300693041e-13
# C_tt_air_out: 1.3685431112559984e-13

# 積分２
f = lambda phi, r: np.sinc(phi/np.pi)/( ((delta_out-delta_in)/(r_out-r_in))*(r-r_in + delta_in + Do*(1-np.cos(phi))))
val, err = integrate.dblquad(f, r_in, r_out, lambda phi: (-1/2)*np.pi, lambda phi: (1/2)*np.pi)

C_tt_air_l = (Do/2.0)*eps_0 * val
print("C_tt_air_l", C_tt_air_l)
# C_tt_air_l 4.721164453460405e-13

C_tt_ins = (np.pi*eps_0*eps_r*L_T)/(2*np.log(Do/Di))
print("C_tt_ins", C_tt_ins)
# C_tt_ins 7.684597167088852e-11

C_tt_in = (C_tt_ins*C_tt_air_in)/(C_tt_ins+C_tt_air_in)
C_tt_out = (C_tt_ins*C_tt_air_out)/(C_tt_ins + C_tt_air_out)
C_tt_l = (C_tt_ins*C_tt_air_l)/(C_tt_ins + C_tt_air_l)

C_tt = C_tt_in + C_tt_out + 2*C_tt_l
print(C_tt)
# 1.6141283310668653e-12

追記

import sympy as sym
sym.var("phi delta_out delta_in r_out r_in r Do")
f = sym.sin(phi)/(phi*( ((delta_out-delta_in)/(r_out-r_in))*(r-r_in + delta_in + Do*(1-sym.cos(phi)))))
display(f)

Comments

[mag0123]

jeanbiego様
ご回答いただき誠にありがとうございます。
修正していただいたコードをこちらでも試したところ、同様の結果を得ました。

コード中の式は論文の式を参考にしていますが、一部は自分が式を設定しているものもあったりしますので一度式の部分を見直してみようと思います。

もう1点気になるところなのですが、積分２の重積分は正しく記述できていますでしょうか。
積分2で実現したい式の画像を補足情報に追加しています。
お忙しいところ恐縮ですが、ご回答いただければ幸いです。

[jeanbiego]

sympyを使ってfを表示してみましたが（回答欄に追記）、r-r_inの括弧位置がおかしいせいで、実現したい式とは異なっていますね。
あと、画像の重積分のところですが、rとφの積分範囲が逆ではないですか？

[mag0123]

jeanbiego様
ご回答いただき誠にありがとうございます。
ご指摘していただいた式を修正したところ、計算結果が1.0206928711387264e-12で、1.124E-12との誤差が約-9.3%とだいぶ近づきました。
見落としていたところをご指摘いただきとても助かりました。

画像の積分範囲も逆になっていました。正しい積分の順序はφ→rの順だと思います。
コード中ではこれが実現できていますでしょうか。
お忙しいところ再々の質問申し訳ございません。

[jeanbiego]

コード中では重積分の範囲等は問題ないように見えます。

[mag0123]

jeanbiego様
ご回答いただき誠にありがとうございます。
重積分の範囲の確認ありがとうございます。

この度は問題の解決にご協力いただきありがとうございました。