问题描述
将 2
到 inf
的跟随函数与 scipy.integrate.quad
集成。
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,radius3))
print('w2=',w2[0])
print('w3=',w3[0])
我使用 n2=22144
和 n3=22145
,但得到两个完全不同的结果。其实这两个结果应该很接近。
w2= 0.4256854383924181
w3= 11774635265351.027
更有趣的是,当将4.11447
(在被积函数中的“return”行)更改为4.1145
或4.1144
时,结果将接近。
#4.1145
w2= 11768644300129.305
w3= 11771237857123.924
#4.1144
w2= 0.42568494248362776
w3= 0.4256854757392548
哪个结果是正确的?
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)