0

将 follow 函数从2toinf与集成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, 2, np.inf, args=(radius1, radius3))
print('w2=', w2[0])
print('w3=', w3[0])

我使用n2=22144and n3=22145，但得到两个完全不同的结果。实际上，这两个结果应该非常接近。

w2= 0.4256854383924181
w3= 11774635265351.027

更有趣的是，当将4.11447（在被积函数中的“返回”行）更改为4.1145or4.1144时，结果将接近。

#4.1145
w2= 11768644300129.305
w3= 11771237857123.924

#4.1144
w2= 0.42568494248362776
w3= 0.4256854757392548

哪个结果是正确的？

4

0 回答 0