Python 牛顿拉夫森小数根

Question

我有这个函数（Newton-Raphson算法）： "digits": 根的期望精度

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 函数有疑问。有什么想法吗？

提前致谢

Answer 1

如果您在每个步骤中打印 x 和 derivative 的值，您将看到：

1.9347284759 0.0683501539811
1.95705673873 0.030810070939
1.97161234728 0.0138147377982
...
1.99997209174 1.42108547152e-08
1.99997909174 7.1054273576e-09
inf 0.0
nan nan

我们可以看到导数为零，除以零为 inf。

如果你使用更大的dx收敛。对于 dx=1e-3

...
1.99991975788 2.71786149142e-06
1.99992026519 2.7172077921e-06
1.99992075954 2.7165576455e-06
Root: 1.999921
Number of loops: N=95

Answer 2

Newton-Raphson 不适用于导数为零的根。正如该图所示，您的问题就是这种情况。

事实上，导数为零的任何根都很难得到高精度，这适用于任何数值求根方法。你的情况更糟，因为一阶和二阶导数在根处都为零。这意味着函数的图形非常接近根部的水平直线。鉴于此，将很难获得比有效数字总数的三分之一更好的精度。大多数计算机会给出双精度，大约 15 位数字，因此超过 5 位有效数字将很困难。您的例程已经给出了，所以不要期望任何数字例程会更好。

当您将导数的精度取为 1e-6（derivative() 调用中的参数）时，您尤其会遇到问题，而当您的 x 值达到相同的精度时，问题就会出现在 x.

您可以使用 Newton-Raphson 的某些变体来避免这个问题和其他问题——请记住，Newton-Raphson 不能保证在一般情况下收敛并且有很多困难。查看 数值食谱 一书中的 rtsafe 作为 N-R (Newton-Raphson) 的安全用法。

如果您想非常接近 N-R，您可以修改代码以在计算 x 的下一个值之前检查导数的值。当导数变为零时停止循环。

如果您真的想在像您这样的情况下获得更高的精度并使用直接 N-R，则需要更好地计算导数。您使用的数值导数显然不是很好。该数值导数可能会给出二阶结果，使用评估点周围两点的对称差异，这在这里还不够好。您可以使用更复杂的衍生程序，它使用更多的点并给出更高阶的结果。但是在您的情况下，该函数非常简单，您可以编写自己的函数，为您的原始函数提供良好的衍生。只需使用基本的微积分规则直至乘积规则。事实上，任何实际的 Newton-Raphson 程序都应该有一个计算导数的函数，除非必要，否则您通常不使用近似的数值导数。有人可能会争辩说，使用数值导数意味着您实际上并不是在做牛顿-拉夫森，您应该使用另一种方法来求根——Numeric Recipes 一书实际上说明了这一点.

修改您的代码以使用定义的导数函数，以及清理您的一些样式，

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

Python 牛顿拉夫森小数根

Python Newton Raphson Decimal Roots

python

precision

decimal