24小时热门版块排行榜    

查看: 2190  |  回复: 4

swuli

木虫 (小有名气)

[求助] scipy.integrate.quad 积分错误原因? 已有1人参与

有以下积分,n2和n3的值非常接近,但积分结果却大相径庭。
n2=22144  积分结果=0.4256854383924181
n3=22145 积分结果=11774635265351.027

用matlab算,后面是准确的。是什么原因导致n2所出现的错误?


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])
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
2楼2021-05-21 16:59:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bcsnow

铁杆木虫 (著名写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与,应助指数 +1
swuli: 金币+10, ★★★★★最佳答案 2021-05-23 08:45:33
没有错误啊:w2= 11772040939067.406  w3= 11774635265351.027
3楼2021-05-22 10:36:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

swuli

木虫 (小有名气)

引用回帖:
3楼: Originally posted by bcsnow at 2021-05-22 10:36:12
没有错误啊:w2= 11772040939067.406  w3= 11774635265351.027

我这里得到的结果却是w2= 0.4256854383924181;w3= 11774635265351.027。奇怪了。
和版本有关吗?我的python 3.6.6 scipy 1.6.0。您有没有做什么别的处理呢?
4楼2021-05-23 08:48:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bcsnow

铁杆木虫 (著名写手)

引用回帖:
4楼: Originally posted by swuli at 2021-05-23 08:48:12
我这里得到的结果却是w2= 0.4256854383924181;w3= 11774635265351.027。奇怪了。
和版本有关吗?我的python 3.6.6 scipy 1.6.0。您有没有做什么别的处理呢?...

没有处理,直接拷贝运行的,我的是3.7.7和1.5.2,按说版本间不应该有这么严重错误。你把n2 n3的值调换一下,看看结果如何。
5楼2021-05-24 00:23:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 swuli 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[硕博家园] 极限拉扯,如何抉择 +7 Isak0320 2026-10-12 10/500 2026-10-12 22:38 by 大飞鱼鱼鱼
[基金申请] 入职后第一次管课题组的专利缴费,对账才发现官费比别人多掏了一倍 +6 13108017953 2026-10-08 6/300 2026-10-12 20:07 by w3w3y3
[有机交流] OpenAI又攻克数学难题了,离AI做有机合成还有多远? +4 wangyl_123 2026-10-10 8/400 2026-10-12 18:09 by wangyl_123
[教师之家] 课题组会帮你通过考核么? +3 好运来~9527 2026-10-10 6/300 2026-10-12 14:55 by 简单化xn
[公派出国] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 22:34 by XQNyDCa1Tk5n
[找工作] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 22:26 by XQNyDCa1Tk5n
[公派出国] 售SCI一区文章,我:8O5.5.1.O5.4,科目全,可伽急 +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 22:26 by XQNyDCa1Tk5n
[考博] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 9/450 2026-10-11 22:15 by XQNyDCa1Tk5n
[考博] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 22:11 by XQNyDCa1Tk5n
[教师之家] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 (EPI+-1)(金币-50) +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 22:06 by XQNyDCa1Tk5n
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 (EPI+-1)(金币-50) +5 qR2AdCWTVtkw 2026-10-09 10/500 2026-10-11 14:25 by VOMexwnDzOFZ
[教师之家] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 (EPI+-1)(金币-50) +4 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 14:05 by VOMexwnDzOFZ
[考研] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +5 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 13:51 by VOMexwnDzOFZ
[硕博家园] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +5 qR2AdCWTVtkw 2026-10-09 8/400 2026-10-11 13:45 by VOMexwnDzOFZ
[找工作] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 10:54 by VOMexwnDzOFZ
[找工作] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +4 qR2AdCWTVtkw 2026-10-09 7/350 2026-10-11 10:34 by VOMexwnDzOFZ
[考博] 招收博士生一名 +3 mandolin 2026-10-10 3/150 2026-10-11 10:21 by 不是山谷!
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 qR2AdCWTVtkw 2026-10-09 4/200 2026-10-11 03:14 by M0zrYVSAnxpE
[有机交流] 各位大佬,路线设计求助 +6 zapen 2026-10-08 8/400 2026-10-10 15:59 by zapen
[有机交流] 我应该选哪条合成路线让实验成功率高一些? +5 jing816y 2026-10-08 9/450 2026-10-10 09:05 by zapen
信息提示
请填处理意见