这是两份代码,一份是用 Python 3 编写的,另一份是用 Wolfram 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()
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 (正确一个)
结果不同。有人知道为什么吗?
最佳答案
从我的测试来看,一阶 Bessell 函数似乎给出了不同的结果。两者最初都评估为 Bessel(f * 0.0188925),但 scipy 版本给出了 0 到 9.4e-3 的范围,其中 wolframalpha(使用 Mathematica 后端)给出 0 到 1.4。我会对此进行更深入的研究。
此外,python 使用标准 C float ,而 Mathematica 使用符号运算。 Sympy试图在 Python 中模仿这种符号操作。
关于python-3.x - 等效代码,不同结果(Python、Mathematica),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/45092918/