Pythonでニュートン法ニューモン - 工作の２S(Seiri/Seiton)

●あらまし

高校のプログラミングの授業に発展問題としてこの「ニュートン法」が出ました
歯抜けのフローチャートがあり選択式で埋めるという問題でした

今思えば情報技術者試験の引用とすぐ分かるものですが
まったく手も足も出ずという感じだったのを覚えています
Python入門のタイミングでリベンジしたら解けたので
まとめることにしたのです

●概要

今回は情報技術者試験に登場したニュートン法
放物線　f(x) = x^2 - N = 0となる解を１つ求めるというものです
f:id:kt-marshal:20160822000055j:plain

●考え方

十分に大きい数を選びx0とする
そのときのf(x)の接線とx軸との交点のx座標をx1とする
x1をx0とし、同様に接線からx1を求める　・・・くりかえし
f:id:kt-marshal:20160822220406j:plain
f:id:kt-marshal:20160822220413j:plain

ひたすら繰り返すとx0 ≒ x1 となる（収束）
プログラムでは欲しい精度を満たす段階で打ち切る
　|x0 - x1| < ε

●Pythonで実装

# coding: utf-8
print("solve \"f(x) = x^2 - N \"")

N = float(input("input N->"))
ep = float(input("input convergence value->"))
x1 = float(input(input start x value->"))
x0 = 0.0
N /= 2

while(abs(x0 - x1) >= ep):
    x0 = x1
    x1 = (x0 / 2) + (N / x0)
    print("x0->{0}, x1->{1}".format(x0, x1))
    print(x1)
    print("")

print(x1)

▼N = 10
f:id:kt-marshal:20160906213556j:plain

●解説

情報技術者試験では接線を求める公式を当てはめるような気がします
今回で行けば与えられた式から「x0を代入するだけで次のx1が求まる式」を導くことになります

f(X)における接線yは
  Y = f'(X)(Xn - X) + f(X)

XをX0, XnをX1と置き換える
  Y = f'(X0) * (X1 - X0) + f(X0)
f'(X) = 2Xより　f'(X) = 2X0
したがって
  y = 2X0(X1 - X0) + f(X0)
    = 2*X0*X1 - 2X0 + f(X0)
よって
  x1 = (X0 / 2) / (N / 2X0)

プログラム内の
　x1 = (x0 / 2) + (N / x0)
が求めた式に相当します
（Nを毎回2で割ることが分かっているためあらかじめ割っています）
f:id:kt-marshal:20160906214019j:plain