我有这个功能（牛顿-拉夫森算法）：
“数字”：所需的根精度
from scipy.misc import derivative
def newtonDigits(function,xstart,digits):
xprev=0
ncalls=0
while abs(xstart-xprev) >= 0.5 * 10**(-digits):
ncalls +=1
print xstart
x = xstart - (function(xstart)/derivative(function,xstart,dx=1e-6))
xprev = xstart
xstart = x
return xstart,ncalls
输入为：
f = lambda x: 14*x*(math.e)**(x-2) - 12*(math.e)**(x-2) - 7 *(x**3) + 20*(x**2) - 26*x + 12
root = newtonDigits(f,1.9,6)
print "Root: {0:.6f}".format(root[0])
print "Number of loops: N=" + str(root[1])
输出为：
1.9
1.9347284759
1.9570567413
1.97161234354
1.98117864941
1.98749755734
1.99168484127
1.99446528727
1.99631400887
1.99754426317
1.9983635866
1.99890940938
1.9992730728
1.99951577542
1.99967709739
1.99978645669
1.99985605669
1.99990268169
1.9999354436
1.9999606936
1.9999726936
1.9999806936
1.9999846936
newtonr.py:55: RuntimeWarning: divide by zero encountered in double_scalars
x = xstart - (function(xstart)/derivative(function,xstart,dx=1e-6))
mainfile.py:6: RuntimeWarning: invalid value encountered in double_scalars
f = lambda x: 14*x*(math.e)**(x-2) - 12*(math.e)**(x-2) - 7 *(x**3) + 20*(x**2) - 26*x + 12
inf
Root: nan
Number of loops: N=24
我尝试使用Decimal，但是浮点数和lambda函数有问题。
有任何想法吗？
提前致谢
最佳答案
牛顿-拉夫森不适用于导数为零的根。如图所示，您的问题就是这种情况。
实际上，导数为零的任何根都很难高精度地获得，并且对于任何数值根查找方法都成立。您的情况更糟，因为一阶和二阶导数的根都为零。这意味着函数的图非常靠近根部的水平直线。鉴于此，将很难获得比有效数字总数的三分之一高的精度。大多数计算机将提供双精度，大约15位数字，因此要超过5位有效数字将非常困难。您的例程已经给出了，所以不要期望任何数字例程都会有更好的结果。
当您将导数的精度带到1e-6（您的derivative()调用中的参数）时，您尤其会遇到问题，而当x值达到x中相同的精度时，就会出现问题。
您可以使用Newton-Raphson的某些变体来避免此问题和其他问题-请记住，Newton-Raphson不能保证在一般情况下收敛，并且有很多困难。请看《数值食谱》一书中的rtsafe作为N-R（牛顿-拉夫森）的安全用法。
如果要非常靠近N-R，可以在计算下一个x值之前修改代码以检查导数的值。当导数变为零时，只需停止循环。
如果您真的想在像您这样的情况下使用更高的精度并使用直线N-R，则需要对导数进行更好的计算。您使用的数字导数显然不是很好。该数值导数可能使用评估点周围两个点的对称差给出二阶结果，在此还不够好。您可以使用更复杂的派生例程，该例程使用更多的点并给出更高阶的结果。但是在您的情况下，函数足够简单，因此您可以编写自己的函数，该函数可以很好地继承原始函数。只需使用基本演算规则，直到乘积规则为止。实际上，应该给任何实际的Newton-Raphson例程一个计算导数的函数，除非必要，否则通常不应使用近似的数字导数。有人可能会争辩说，使用数值导数意味着您实际上并没有在做牛顿-拉弗森，而应该使用另一种方法来找到根-《数字食谱》一书确实指出了这一点。
修改代码以使用定义的派生函数，以及清理某些样式，可以
import math
def newtonDigits(function, dfunction, xstart, digits):
xprev=0
ncalls=0
while abs(xstart-xprev) >= 0.5 * 10**(-digits):
ncalls +=1
print(xstart)
x = xstart - function(xstart) / dfunction(xstart)
xprev = xstart
xstart = x
return xstart, ncalls
f = lambda x: 14*x*(math.e)**(x-2) - 12*(math.e)**(x-2) - 7*x**3 + 20*x**2 - 26*x + 12
df = lambda x: 14*x*(math.e)**(x-2) + 2*(math.e)**(x-2) - 21*x**2 + 40*x - 26
root = newtonDigits(f, df, 1.9, 6)
print("Root: {0:.6f}".format(root[0]))
print("Number of loops: N=" + str(root[1]))
哪个打印
1.9
1.9347284749567717
1.9570567399626235
1.9716123428786936
1.9811786535860885
1.987497587688875
1.991684850749086
1.9944652840851569
1.9963140403238206
1.9975443983678638
1.9983636880398847
1.998909460162945
1.9992731216554642
1.9995154801019759
1.9996770218573676
1.999784687490283
1.9998564592822783
1.9999043185996854
1.9999363386081155
1.9999575985093117
1.9999719257797395
1.999982068293325
1.9999875928805002
2.00000490230375
Root: 2.000005
Number of loops: N=24