24小时热门版块排行榜    

查看: 1016  |  回复: 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 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[论文投稿] 《控制理论与应用》期刊收版面费吗? +5 ygj2015 2024-05-29 5/250 2024-06-03 09:19 by lmpmsm
[论文投稿] Scientific reports 投稿 +4 chencome12 2024-06-02 6/300 2024-06-03 09:13 by chencome12
[基金申请] 前些天开会有个人见到人就搞关系,一查此人全是MDPI/Hindawi论文,鄙视! +9 zju2000 2024-06-02 9/450 2024-06-03 09:03 by xiaoding
[基金申请] 化学口B0109(高分子合成),拿青年基金一般需要怎样的文章水平? +14 salmon95 2024-05-30 24/1200 2024-06-03 08:31 by LNP@mRNA
[论文投稿] 选期刊 5+3 jfdhj 2024-05-29 5/250 2024-06-03 07:47 by 黑大环境队长
[基金申请] 2024杰青和万人领军什么时候会评 +6 墨香琴韵 2024-06-02 6/300 2024-06-03 07:24 by llhljsy
[硕博家园] 又想换工作 +15 brightmj 2024-05-27 26/1300 2024-06-03 07:01 by zyqchem
[论文投稿] 求助大神,Fe和Al离子对MOF都有淬灭,当两种离子共存时,怎么区分两种离子? 10+3 maoxiao 2024-06-02 3/150 2024-06-03 01:43 by nono2009
[论文投稿] 希望嫡长子顺利出场! +5 C刊霸王 2024-05-31 6/300 2024-06-02 19:26 by 子非愚
[基金申请] 入职高校3年发表10+SCI,尽人事听天命 +33 kaoyan250 2024-05-27 48/2400 2024-06-02 19:17 by kaoyan250
[考博] 导师不让硕转博,让我去国外读博,能理解吗? +12 萧山幽谷 2024-05-29 20/1000 2024-06-02 12:01 by yuan0806
[考博] 申请2024或2025年博士研究生 +5 嘟噜嘟1 2024-05-29 12/600 2024-06-01 22:36 by 嘟噜嘟1
[考博] 24年博士招生 +8 abinit432 2024-05-27 10/500 2024-06-01 17:38 by czp97
[考博] 24or25材料专业申博 +4 农夫三拳有点痛 2024-05-30 11/550 2024-06-01 14:45 by Napoleonsky
[基金申请] 博后特别资助状态变化 +24 随梦而飞2017 2024-05-30 35/1750 2024-06-01 10:10 by 青岛阳仔
[硕博家园] 哈工大硕博招生!博士每月入学! +4 nailooo 2024-05-30 5/250 2024-06-01 06:47 by anevay
[高分子] MMA预聚体光固化发雾问题求助 +3 惠亚金总 2024-05-29 10/500 2024-05-31 14:59 by 惠亚金总
[论文投稿] 有没有老师需要发表论文 +4 金老师论文助理- 2024-05-29 4/200 2024-05-29 16:51 by liuyupu132
[基金申请] E05青基有几个评审 +4 KYXY123 2024-05-28 4/200 2024-05-28 19:25 by popt2t
[基金申请] E10开始送了,希望有好运 +5 sail 2024-05-27 5/250 2024-05-28 18:36 by 芝小芝
信息提示
请填处理意见