解决精度问题 (python)

Question

我遇到了 fsolve 的精度问题。

import numpy as np
from scipy.optimize import fsolve


ens=np.arange(0,50,1)


def f(x):
     return x*(np.sin(x)/np.cos(x))-1000

s=[]        

roots=fsolve(f,ens)
roots=np.around(roots, decimals=3 , out=None)
a = roots[roots >= 0]
g = np.unique(a)
g=g[:5]
s.append(g)
print(s)
            

result :

[array([10.842, 11.006, 15.165, 21.116, 22.382])]

结果应该是：[1.569,4.708,7.846,10.985,14.123]

我的代码遗漏了前三个解决方案，其他的也不准确。您知道我怎样才能提高结果的精确度吗？

Answer 1

您可以使用 scipy.optimize.newton() with first and second derivative of the function f to get a more accurate result. You can do it by hand, or use derivative-calculator.net or wolframalpha. Pass both the first derivative fprime and the second derivative, fprime2 to newton(). If you do this, Halley's method will be used instead of a simple Newton-Raphson, which is more accurate (see the description below fprime2 in newton() 文档）。

def f(x):
     return x*(np.sin(x)/np.cos(x))-1000

def fprime(x):
    # first derivative of f
    # do this by hand or use derivative-calculator.net or wolframalpha
    return np.tan(x)+ x*(1/np.cos(x))**2

def fprime2(x):
    # second derivative of f
    # do this by hand or use derivative-calculator.net or wolframalpha
    return 2*1/(np.cos(x)**2)*(x*np.tan(x)+1)

res = newton(f,x0=[1,4,7,10,14],fprime=fprime, fprime2=fprime2)
res = np.around(res, decimals=3)

res 将是：

array([ 1.569,  4.708,  7.846, 10.985, 14.123])

上面newton()中的x0参数是所谓的根所在的初始猜测的列表。由于您的函数有无限多个根（见下图），传递其中一些有助于获得您真正关心的根。

这就是 f 的样子（好吧，f/1000 让特征可见）：

解决精度问题 (python)

Fsolve precision issue (python)

python

numpy

precision

solver