2011-08-16 71 views
0

我有解決牛頓方法的代碼,在上一次迭代中,輸出值出錯了。這些值不應該是負值,因爲我在紙上手動檢查了它。據我知道的代碼是正確的,但我不能在0和1之間。這裏是代碼弄清楚爲什麼會顯示負值,也是最終的U值最好應爲正值:牛頓方法的意外輸出

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() 

對於最後一次迭代的輸出是下面,

其中K6 = -0.02,K5 = 0.03,K4 = -0.03和K0 = 0.0

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 

的初始猜測值是0.5,所以如果它是取代的在多項式中,那麼輸出是-0.005。

然後再次在微分多項式中替換相同的初始值。結果是-0.015。

現在這些值在牛頓方程中被代入,那麼答案應該是0.166666667。但實際答案是一個負值。

謝謝。

+1

** 1 **:你能解釋'poly_coeff'嗎?該算法來自哪裏?我所知道的所有其餘數學都是正確的,並且我得到了你發佈的答案。 ** 2 **:你說你已經在紙上手動檢查過了 - 你可以發佈這些計算,所以我們可以看到預期和實際之間有什麼區別? – Nate

+0

@注意poly_coeff()計算註釋中給出的多項式的係數。然後用牛頓法計算該多項式的根。我將編輯顯示我的手動計算的問題。謝謝 – zingy

回答

0

啊,我明白了。

正如你說的,

float(poly_substitute(poly, x)) 

評估爲-0.015。然後,

float(poly_substitute(poly_diff, x)) 

評估爲-0.01。因此,替代這些值和x

x_n = 0.5 - ((-0.015)/(-0.01)) 
x_n = 0.5 - 0.6666666... 
x_n = -0.166666... 

你手動數學是什麼過錯,而不是代碼。

+0

是的,我發現代碼很好。另外它的python浮點問題給了意想不到的結果。謝謝 – zingy

0

給出的多項式在x = 0時有一個單一的解。代碼工作正常。