ニュートン法

Pythonでニュートン法をやってみた。ニュートン法って何？って方はググってください。

def newton(x, num):
    #nextx = (2 * (x*x*x) + x*x - 3) / (3 * (x*x) + 2*x - 5)    #x - (g / dg)をまとめた場合
    g = f(x)
    dg = df(x)
    nextx = x - (g / dg)
    num.append(nextx)
    if x != 0:
        check = (nextx - x) / x
        if (check < 0.00000001 and check > -0.00000001):    #誤差の範囲内に収まったら計算終了
            return
    newton(nextx, num)

def f(x):
    return x*x*x + x*x - 5*x + 3    #今回はx^3+x^2-5x+3=0を解く

def df(x):
    return 3*x*x + 2*x - 5    #左辺を1回微分した式

def main():
    x0 = 0.0    #最初のxの値
    num = [x0]
    newton(x0, num)
    print num

if __name__ == '__main__':
    main()

結果は以下のようになる。

>>> 
[0.0, 0.6, 0.8117647058823529, 0.9082650781831724, 0.9546772328747365, 0.9774692207697667, 0.9887666079847501, 0.9943912241748202, 0.9971975823794107, 0.9985992825526064, 0.9992997639662882, 0.9996499126368127, 0.9998249639795537, 0.9999124839046458, 0.9999562424305699, 0.9999781213351223, 0.9999890606966689, 0.9999945303568718, 0.9999972651704089, 0.9999986325565382, 0.9999993162552172, 0.9999996580523504, 0.9999998288321135, 0.9999999147737951, 0.9999999577621707, 0.9999999801044548, 0.999999988474844]

最初のxの値を5.0にした場合は以下のようになる。

>>> 
[5.0, 3.4, 2.389473684210526, 1.774066654525587, 1.4160571394231234, 1.2173873680856377, 1.1114246354069843, 1.056457160596543, 1.0284236632976032, 1.0142617931551174, 1.007143541372518, 1.0035749515517178, 1.0017882734738675, 1.000894336473054, 1.0004472182097137, 1.000223621603005, 1.0001118139264726, 1.0000559077440478, 1.0000279540660224, 1.000013977080193, 1.0000069885504108, 1.0000034942696048, 1.000001747124112, 1.0000008735613346, 1.0000004367476691, 1.0000002183878496, 1.0000001093420305, 1.0000000545122083, 1.0000000270174305, 1.0000000105802935, 1.0000000105802935]