Python牛顿·拉夫森;计算精度

Python newton raphson; precision in calculations

我在python中写了这个函数:

def f2(x):
    return (5.0*x + log1p(x) - 10000.0)

def dfdx2(x):
    return (5.0-(1.0/x))

def newtonRaphson2(f, dfdx, x, tol):

    x0 = x

    for i in range(1, 2000):
        if f(x) == 0.0:
            return x

    if dfdx(x) == 0.0:
        print(dfdx(x))
        break

    x = x - (f(x) / dfdx(x))
    #print(x)

    Er = abs(x0-x)/abs(x0)
    if Er <= tol:
        return x
    print(Er)
    x0 = x
    return x

然后我这样执行:

task2 = newtonRaphson2(f2, dfdx2, 1, 0.000001);
print(task2)

对于输出检查 Er 在它之前输出的最终精度为 4.245665128208564e-05 returns x.

X 在 1998.479871524306 处返回,这是一个很好的估计,但我最好至少将其降低到 1.0e-06。将 tol 变量更改为 1.0e-08 似乎没有任何作用。

我猜也许将每个变量都变成双精度是个更好的主意,但我仍然不知道为什么我的代码停在它该做的地方。我对 python 也不是那么稳定,这就是我问的原因。我已经编写了其中一个有效的方法,但它用于一个简单得多的方程式。

正如@alex-dubrovsky提到的,意味着收敛的计算应该用条件循环来实现,即:

while True:
    if f(x) == 0.0:
        return x

    if dfdx(x) == 0.0:
        print(dfdx(x))
        break

    x = x - (f(x) / dfdx(x))
    #print(x)

    Er = abs(x0-x)/abs(x0)
    if Er <= tol:
        return x
    print(Er)
    x0 = x

使用这种方法,您总是面临无限循环的风险,但这或多或少没问题,正如算法所暗示的那样"running until converged"

你的代码工作正常,一旦你正确缩进并添加 from math import log1p 。只需将 print(Er) 行放在计算后的右边即可查看其最终值。 Er 达到 ~10^-9。这对我有用:

from math import log1p

def f2(x):
    return (5.0*x + log1p(x) - 10000.0)

def dfdx2(x):
    return (5.0-(1.0/x))

def newtonRaphson2(f, dfdx, x, tol):

    x0 = x
    for i in range(1, 2000):
        if f(x) == 0.0:
            return x

        if dfdx(x) == 0.0:
            print(dfdx(x))
            break
        x = x - (f(x) / dfdx(x))
        #print(x)

        Er = abs(x0-x)/abs(x0)
        print('Er = {}'.format(Er))
        if Er <= tol:
            return x

        x0 = x
    return x

x = newtonRaphson2(f2, dfdx2, 1, 0.0000001)
print 'X = {}'.format(x)

输出为:

Er = 2498.5767132
Er = 0.200506616666
Er = 4.24566512821e-05
Er = 8.49642413214e-09
X = 1998.47987152

考虑在此处使用 while。 Newton-Raphson 算法通常收敛得非常快,因此在大多数情况下 运行 不需要多次迭代。

这给出了相同的结果:

    from math import log1p

def f2(x):
    return (5.0*x + log1p(x) - 10000.0)

def dfdx2(x):
    return (5.0-(1.0/x))

def newtonRaphson2(f, dfdx, x, tol):

    x0 = x
    Er = 1
    while Er >= tol:
        if f(x) == 0.0:
            return x

        if dfdx(x) == 0.0:
            print(dfdx(x))
            break
        x = x - (f(x) / dfdx(x))
        #print(x)

        Er = abs(x0-x)/abs(x0)
        print('Er = {}'.format(Er))

        x0 = x
    return x

x = newtonRaphson2(f2, dfdx2, 1, 0.0000001)
print 'X = {}'.format(x)