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]
また、グラフは以下のようになる。
それっぽい値に収束していることが分かる。