为啥 sympy 在整合部分功率时给出不同/错误的答案？

Posted 2023-03-24

【中文标题】为啥 sympy 在整合部分功率时给出不同/错误的答案？

【英文标题】：Why sympy gives different/wrong answer when integrate fraction of power?为什么 sympy 在整合部分功率时给出不同/错误的答案？

【发布时间】：2021-12-19 23:32:02

【问题描述】：

我正在运行一个定积分，但 sympy 给了我一个不同的答案。 我不知道出了什么问题？ 任何人都可以帮忙吗？谢谢

The equation integrated by Wolfram, the result is as expected 1.097

import sympy as sp
y=sp.Symbol('y')
fy=(y-2.4)**2*sp.sqrt(y)*0.1875
print(sp.integrate(fy,(y,0,4)))

这给了我 0.485255751108899

【问题讨论】：

对于它的价值，scipy.integrate.quad 同意 Wolfram。

@TimRoberts, quad(f,2.4,4) 是 0.485，sp.integrate 也是如此。出于某种原因，sympy 错过了第一个驼峰，在 (0,2.4) 范围内返回 0。更具体地说，sympy 首先上升，(0, .5) 错误。

integrate(fy,(y,0.001,4)) 是 1.0971。出于某种原因，从 0 开始集成会给它带来问题。也许是因为绘图线在 0 处是垂直的。数值积分不会受到影响。 @蒂姆罗伯茨

【参考方案1】：

这看起来像是一个与浮点数有关的错误，在集成中可能会出现问题：

In [15]: from sympy import *

In [16]: y = symbols('y')

In [17]: fy=(y-2.4)**2*sp.sqrt(y)*0.1875

In [18]: fy
Out[18]: 
                                 2
1.08⋅√y⋅(0.416666666666667⋅y - 1) 

In [19]: integ = Integral(fy, (y, 0, 4))

In [20]: integ
Out[20]: 
4                                      
⌠                                      
⎮                                  2   
⎮ 1.08⋅√y⋅(0.416666666666667⋅y - 1)  dy
⌡                                      
0   

In [21]: integ.doit() # symbolic integration with floats
Out[21]: 0.485255751108899

In [22]: integ.evalf() # numeric integration
Out[22]: 1.09714285714286

In [23]: integ_rational = Integral(nsimplify(fy), (y, 0, 4))

In [24]: integ_rational # same integral with Rational instead of Float
Out[24]: 
4                    
⌠                    
⎮                2   
⎮       ⎛5⋅y    ⎞    
⎮ 27⋅√y⋅⎜─── - 1⎟    
⎮       ⎝ 12    ⎠    
⎮ ──────────────── dy
⎮        25          
⌡                    
0                    

In [25]: integ_rational.doit() # symbolic integration with rationals
Out[25]: 
192
───
175

In [26]: _.n() # numerically evaluate the last result
Out[26]: 1.09714285714286

首先不要使用带有 sympy 的浮点数，除非你有充分的理由，例如使用Rational('0.1') 或S(1)/10 等而不是0.1。 nsimplify 函数可用于将 float 转换为 Rational（它会尝试猜测您的真正意思是什么）。

其次，这是一个错误，因此应将其报告给 github 而不是 SO： https://github.com/sympy/sympy/issues

【讨论】：

您好，奥斯卡，感谢您的回复。我会向github报告。

github.com/sympy/sympy/issues/22435

使用nsimplify(2.4) 足以让这个工作正常。