牛顿法的意外输出

使用者使用牛顿法编写了这段代码，但最后一步的输出值是错误的。他手动计算了结果，确信输出值不会是负数。代码如下：

import copy
import math

tlist = [0.0, 0.07, 0.13, 0.15, 0.2, 0.22] # list of start time for the phonemes

w = 1.0

def time() :
    t_u = 0.0
    for i in range(len(tlist)- 1) :
        t_u = t_u + 0.04 # regular time interval
        print t_u
        print tlist[i], ' ' , tlist[i+1], ' ', tlist[i -1]
        if t_u >= tlist[i] and t_u <= tlist[i + 1] :
            poly = poly_coeff(tlist[i], tlist[i + 1], t_u)
            Newton(poly) 
        else :
            poly = poly_coeff(tlist[i - 1], tlist[i], t_u)
            Newton(poly) 

def poly_coeff(start_time, end_time, t_u) :
    """The equation is k6 * u^3 + k5 * u^2 + k4 * u + k0 = 0. Computing the coefficients for this polynomial."""
    """Substituting the required values we get the coefficients."""
    t0 = start_time
    t3 = end_time
    t1 = t2 = (t0 + t3) / 2
    w0 = w1 = w2 = w3 = w
    k0 = w0 * (t_u - t0)
    k1 = w1 * (t_u - t1)
    k2 = w2 * (t_u - t2)
    k3 = w3 * (t_u - t3)
    k4 = 3 * (k1 - k0)
    k5 = 3 * (k2 - 2 * k1 + k0)
    k6 = k3 - 3 * k2 + 3 * k1 -k0 

    print k0, k1, k2, k3, k4, k5, k6
    return [[k6,3], [k5,2], [k4,1], [k0,0]]

def poly_differentiate(poly):
    """ Differentiate polynomial. """
    newlist = copy.deepcopy(poly)

    for term in newlist:
        term[0] *= term[1]
        term[1] -= 1

    return newlist

def poly_substitute(poly, x):
    """ Apply value to polynomial. """
    sum = 0.0 

    for term in poly:
        sum += term[0] * (x ** term[1])
    return sum

def Newton(poly):
    """ Returns a root of the polynomial"""
    x = 0.5  # initial guess value 
    epsilon = 0.000000000001
    poly_diff = poly_differentiate(poly) 

    while True:
        x_n = x - (float(poly_substitute(poly, x)) / float(poly_substitute(poly_diff, x)))

        if abs(x_n - x) < epsilon :
            break
        x = x_n
        print x_n
    print "u: ", x_n
    return x_n

if __name__ == "__main__" :
    time()

输出结果如下：

0.2
0.2   0.22   0.15
0.0 -0.01 -0.01 -0.02 -0.03 0.03 -0.02
-0.166666666667
-0.0244444444444
-0.000587577730193
-3.45112269878e-07
-1.19102451449e-13
u: -1.42121180685e-26

2、解决方案

问题在于牛顿法的迭代中，迭代值不应该被负值取代。只需添加一个条件检查即可解决此问题，以确保在迭代期间 x_n 不为负值，代码中的检查部分如下：

def Newton(poly):
    """ Returns a root of the polynomial"""
    x = 0.5  # initial guess value 
    epsilon = 0.000000000001
    poly_diff = poly_differentiate(poly) 

    while True:
        x_n = x - (float(poly_substitute(poly, x)) / float(poly_substitute(poly_diff, x)))

        if abs(x_n - x) < epsilon :
            break
        # Ensure x_n is not negative
        if x_n < 0:
            x_n = 0
        x = x_n
        print x_n
    print "u: ", x_n
    return x_n

修正后的代码中的if语句确保在迭代期间 x_n 不为负值，确保Newton方法返回的根是正值。

使用修正后的代码，该使用者得到了正确的结果。