2変数のニュートン法

今度は2変数についてのニュートン法をやってみた。今度は2変数のテイラー展開を利用する。詳しくはググってください。
今回用いた式は某教科書に例として載っていたものをそのまま利用した。

import math
import sys

chk = 0.00000001    #誤差の範囲
count = 0    #計算回数を数える変数

def newton2(x, y, numx, numy):
    global chk, count
    count += 1
    det = fx(x, y) * gy(x, y) - fy(x, y) * gx(x, y)
    nextx = x + ((g(x, y) * fy(x, y) - f(x, y) * gy(x, y)) / det)
    nexty = y + ((g(x, y) * fx(x, y) - f(x, y) * gx(x, y)) / (-1 * det))
    numx.append(nextx)
    numy.append(nexty)
    if x != 0 and y != 0:
        if math.fabs((nextx-x)/x) < chk and math.fabs((nexty-y)/y) < chk:
            return    #x,yが誤差の範囲内になったら終了
        if count > 500:    #500回計算しても収束しなかったら終了
            print 'Error.'
            sys.exit()
    newton2(nextx, nexty, numx, numy)

def f(x, y):
    return x*x + y*y - 1    #x^2+y^2-1=0

def fx(x, y):
    return 2*x    #xについて微分したもの

def fy(x, y):
    return 2*y    #yについて微分したもの

def g(x, y):
    return y - math.sin(x)    #y-sin(x)=0

def gx(x, y):
    return -1 * math.cos(x)    #xについて微分したもの

def gy(x, y):
    return 1    #yについて微分したもの

def main():
    x0 = 1.0    #xの初期値
    y0 = 0.0    #yの初期値
    numx = [x0]
    numy = [y0]
    newton2(x0, y0, numx, numy)
    print numx
    print numy

if __name__ == '__main__':
    main()

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

>>> 
[1.0, 1.0, 0.7566170403394795, 0.7396666657539136, 0.7390853202729097, 0.7390851332151941, 0.7390851332151607]    #numx
[0.0, 0.8414709848078965, 0.7099706104943049, 0.6741397407072063, 0.6736122813110472, 0.6736120291832514, 0.6736120291832148]   #numy

初期値をx=1.0,y=-10.0にすると以下のようになる。

>>> 
[1.0, -12.266955241786267, -9.380381787923056, -6.2671041285218365, -4.769236116058618, -2.3695270054714057, -2.0955558507647027, -0.8394016187699775, -0.9223324082707947, -0.7494421043332685, -0.7392511255062533, -0.7390851493752312, -0.7390851332151609, -0.7390851332151606]
[-10.0, -6.326695524178626, 3.0531088290363426, -3.1545915867907217, 1.5137548246336845, 1.1347277640254845, -0.8939091365913807, -1.4947838723499602, -0.7996337509916909, -0.692593155178541, -0.6737699478467745, -0.6736120504063668, -0.6736120291832152, -0.6736120291832148]

また、グラフは以下のようになる。

それっぽい値に収束していることが分かる。