等效代碼，不同的結果（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）

回答

相關問題