scipy.integrate.quad 積分錯(cuò)誤原因?
有以下積分,n2和n3的值非常接近,但積分結(jié)果卻大相徑庭。
n2=22144 積分結(jié)果=0.4256854383924181
n3=22145 積分結(jié)果=11774635265351.027
用matlab算,后面是準(zhǔn)確的。是什么原因?qū)е耼2所出現(xiàn)的錯(cuò)誤?
from math import sqrt, log, exp
import numpy as np
from scipy.integrate import quad
def integrand(x, r1, r2):
a = r1 + r2
b = r1 * r2
R = a * x / 2
R1 = R ** 2 - a ** 2
R2 = R ** 2 - (r1 - r2) ** 2
vdw = -7e-21 * (2 * b / R1 + 2 * b / R2 + log(R1 / R2))
edl = 7.8e-12 * b / a * log(1 + exp(-328774227 * (R - a)))
force = vdw + edl
return exp(force / 4.11447e-21) / x ** 2
n1 = 5010
n2 = 22144
n3 = 22145
radius1 = 1.05 * 10 ** -8 * n1 ** (1 / 1.8)
radius2 = 1.05 * 10 ** -8 * n2 ** (1 / 1.8)
radius3 = 1.05 * 10 ** -8 * n3 ** (1 / 1.8)
w2 = quad(integrand, 2, np.inf, args=(radius1, radius2))
w3 = quad(integrand, 2, np.inf, args=(radius1, radius3))
print('w2=', w2[0])
print('w3=', w3[0])
返回小木蟲查看更多
京公網(wǎng)安備 11010802022153號(hào)
沒有錯(cuò)誤啊:w2= 11772040939067.406 w3= 11774635265351.027
我這里得到的結(jié)果卻是w2= 0.4256854383924181;w3= 11774635265351.027。奇怪了。
和版本有關(guān)嗎?我的python 3.6.6 scipy 1.6.0。您有沒有做什么別的處理呢?
沒有處理,直接拷貝運(yùn)行的,我的是3.7.7和1.5.2,按說版本間不應(yīng)該有這么嚴(yán)重錯(cuò)誤。你把n2 n3的值調(diào)換一下,看看結(jié)果如何
,