Equivalent codes, different results (Python, Mathematica)

Question

These are two codes, one written with Python 3, and the other one written with Wolfram Mathematica. The codes are equivalent, and therefore the results (plots) should be the same. But the codes give different plots. Here are the codes.

The 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()

the python result

The 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}]

the Mathematica result (correct one)

The results are different. Does someone know why?

Maybe. However, the shift of the minima of the function is more than 0.4. I would not expect this from different float handling. — bogoliuber
– bogoliuber, Commented Jul 14, 2017 at 0:56
A good way to find the source of problem would be, take a sub-expression, check their values individually. Do this check recursively to reduce the size of problem. — Kh40tiK
– Kh40tiK, Commented Jul 14, 2017 at 1:15
Also, I think MMA involve symbolic execution. What happens if you use Compile[] in MMA? — Kh40tiK
– Kh40tiK, Commented Jul 14, 2017 at 1:17
I forgot to mention, the Mathematica result is the correct one. So the problem is in the python code. — bogoliuber
– bogoliuber, Commented Jul 14, 2017 at 1:44

Gavin Achtemeier · Accepted Answer · 2017-07-14 02:12:33Z

1

From my tests it looks like the first order Bessell functions give different results. Both evaluate to Bessel(f * 0.0188925) initially, but the scipy version gives me a range from 0 to 9.4e-3 where wolframalpha (which uses a Mathematica backend) gives 0 to 1.4. I would dig a little deeper into this.

Additionally python uses standard C floating point numbers while Mathematica uses symbolic operations. Sympy tries to mimic such symbolic operations in python.

answered Jul 14, 2017 at 2:12

Gavin Achtemeier

3252 silver badges12 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

bogoliuber Over a year ago

Thanks, will look deeper into it.

Collectives™ on Stack Overflow

Equivalent codes, different results (Python, Mathematica)

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related