等效代码，不同的结果（Python，Mathematica）

这是两个代码，一个用Python 3编写，另一个用Wolfram Mathematica编写。代码是等效的，因此结果（图）应该是相同的。但代码给出了不同的情节。这里是代码。等效代码，不同的结果（Python，Mathematica）

的Python代码：

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.special import k0, k1, i0, i1 

k=100.0 
x = 0.0103406 
B = 80.0 

def fdens(f): 
    return (1/2*(1-f**2)**2+f **4/2 
      +1/2*B*k*x**2*f**2*(1-f**2)*np.log(1+2/(B*k*x**2)) 
      +(B*f**2*(1+B*k*x**2))/((k*(2+B*k*x**2))**2) 
      -f**4/(2+B*k*x**2) 
      +(B*f)/(k*x)* 
      (k0(f*x)*i1(f *np.sqrt(2/(k*B)+x**2)) 
      +i0(f*x)*k1(f *np.sqrt(2/(k*B)+x**2)))/ 
      (k1(f*x)*i1(f *np.sqrt(2/(k*B)+x**2)) 
      -i1(f*x)*k1(f *np.sqrt(2/(k*B)+x**2))) 
      ) 

plt.figure(figsize=(10, 8), dpi=70) 
X = np.linspace(0, 1, 100, endpoint=True) 
C = fdens(X) 
plt.plot(X, C, color="blue", linewidth=2.0, linestyle="-") 
plt.show()

the python result

的数学代码：

k=100.;B=80.; 
x=0.0103406; 
func[f_]:=1/2*(1-f^2)^2+1/2*B*k*x^2*f^2*(1-f^2)*Log[1+2/(B*k*x^2)]+f^4/2-f^4/(2+B*k*x^2)+B*f^2*(1+B*k*x^2)/(k*(2+B*k*x^2)^2)+(B*f)/(k*x)*(BesselI[1, (f*Sqrt[2/(B*k) + x^2])]*BesselK[0, f*x] + BesselI[0, f*x]*BesselK[1, (f*Sqrt[2/(B*k) + x^2])])/(BesselI[1, (f*Sqrt[2/(B*k) + x^2])]*BesselK[1,f*x] - BesselI[1,f*x]*BesselK[1, (f*Sqrt[2/(B*k) + x^2])]); 

Plot[func[f],{f,0,1}]

the Mathematica result （正确的）

结果是不同的。有人知道为什么吗？

来源

2017-07-14 bogoliuber

它们处理浮点的方式不同吗？ – T4rk1n

也许吧。但是，函数最小值的偏移大于0.4。我不希望这从不同的浮动处理。 – bogoliuber

寻找问题根源的一个好方法是，取一个子表达式，单独检查它们的值。以递归方式检查以减少问题的大小。 – Kh40tiK

从我的测试中，它看起来像一阶贝塞尔函数给出不同的结果。最初都评价为贝塞尔（f * 0.0188925），但scipy版本给我的范围从0到9.4e-3，其中wolframalpha（使用Mathematica后端）给出0到1.4。我会深入研究一下。

此外，python使用标准C浮点数，而Mathematica使用符号操作。 Sympy试图模仿Python中的这种符号操作。

来源

2017-07-14 02:12:33

谢谢，会更深入地看。 – bogoliuber

等效代码，不同的结果（Python，Mathematica）

回答

相关问题