等価コード、異なる結果（Python、Mathematica）

Question

これはPython 3で書かれたコードと、Wolfram Mathematicaで書かれたコードの2つのコードです。コードは同等であるため、結果（プロット）は同じでなければなりません。しかし、コードは異なるプロットを与える。ここにコードがあります。

ザ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

Mathematicaコード：

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 （正しい）

結果は異なります。なぜ誰かが知っていますか？

Answer 1

私のテストからは、最初のBessell関数が異なる結果をもたらすように見えます。両方とも最初はBessel（f * 0.0188925）と評価されますが、scipyバージョンでは0〜9.4e-3の範囲が与えられ、wolframalpha（Mathematicaバックエンドを使用）では0〜1.4が与えられます。私はこれを少し深く掘り下げます。

さらに、Pythonは標準C浮動小数点数を使用し、Mathematicaは記号演算を使用します。 Sympyは、このようなシンボリック操作をPythonで模倣しようとしています。