Equivalent codes, different results (Python, Mathematica)
These are two codes, one written with Python 3 and the other with Wolfram Mathematica. The codes are equivalent, so the results (graphs) should be the same. But the codes give different graphs. Here are the codes.
Python code:
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()
Mathematica code:
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}]
Mathematica result (fix)
The results are different. Does anyone know why?
source to share
From my tests, it looks like the first order Bessell functions give different results. Both evaluate Bessel first (f * 0.0188925), but the scipy version gives me a range of 0 to 9.4e-3, where wolframalpha (which Mathematica backend uses) gives 0 to 1.4. I would go deeper into this.
Also, python uses standard C floating point numbers, while Mathematica uses symbolic operations. Sympy tries to mimic such symbolic operations in python.
source to share