任务是解决qubic方程 . 这个方程的系数是依赖于x的函数，所以我在每个点求解方程并得到一个根数组 . 问题在于有些点，其中由sympy给出的所有三个根都是虚构的，但是具有实系数的qubic方程是不可能的 . 此外，不同的sympy版本给出了不同的答案 .
import numpy as np
import pylab
import sympy as sp
from sympy import symbols, solve, Eq
def solFinder(xr, k1r, k2r):
E0r=7446.0
C0r=5E-5
E1r=12500
E2r=17500
L = symbols("L")
solution=[]
for i in range(0, len(xr)):
sol = solve(k1r*k2r*L**3 + (k1r + 2*k1r*k2r*xr[i]-k1r*k2r*C0r)*L**2 + (1 + k1r*xr[i] - k1r*C0r)*L - C0r,L)
solution.append(sol[0].evalf())
return np.array(solution)
x_exp = np.array([0,1E-5, 2E-5, 3E-5, 4E-5, 5E-5, 6E-5, 7E-5])
y_exp = np.array([0,0.08958, 0.14816, 0.19588, 0.21315, 0.22352, 0.2322, 0.24061])
solFinder(x_exp, 10000, 10000)
由sympy 0.7.2-1给出的答案：
array([5.00000000000000e-5, 4.48432635513462e-5, 4.06272178309238e-5,
3.71359664093131e-5, 3.42067801434226e-5, 3.17182646506781e-5,
2.95796710625189e-5, -0.000108861294879183 + 7.86451048332175e-5*I], dtype=object
在最后一点，所有三个根都是虚构的 . 这是最后一点的答案（在evalf（）之前）：
[-0.000108861294879183 + 7.86451048332175e-5*I, -6.33333333333333e-5 - 1.11111111111111e-11*(4092018.02962894 - 7087583.13280521*I)/(-0.5 - 0.866025403784439*I) - (-0.5 - 0.866025403784439*I)*(6.10945499726134e-8 + 1.05818864618122e-7*I), -6.33333333333333e-5 - (-0.5 + 0.866025403784439*I)*(6.10945499726134e-8 + 1.05818864618122e-7*I) - 1.11111111111111e-11*(4092018.02962894 - 7087583.13280521*I)/(-0.5 + 0.866025403784439*I)]
sympy 0.7.1.rc1-3给出的答案是完全不同的：
array([5.00000000000000e-5, -5.74216317756734e-5 + 8.86154779863473e-5*I,
-6.53136089154621e-5 - 8.96724758434974e-5*I,
-7.35679832046565e-5 + 8.97317566604903e-5*I,
-8.21033900717115e-5 - 8.87469355254471e-5*I, 3.17182646506772e-5,
-9.9789835531258e-5 + 8.33396049033042e-5*I, 2.77225886328630e-5], dtype=object)
并且在同情0.7.1.rc1-3中，没有任何一点是所有的根都是虚构的，尽管真正的根是浮动的 .