🐺

平方根をニュートン法により求めてみる

に公開

アドベントカレンダー(自称)vol.16

皆さん、こんにちは。yumyum116 です。SWE転職を目指す者です。
アドベントカレンダーに便乗して、1日1記事投稿に挑戦するとともに、執筆活動の習慣化にも挑戦してみます。

本記事は、頭休めになるような、それでいて、学びにもなるtipsを取り扱います。
今回のテーマは ニュートン法 です。
コーヒー片手にお付き合いいただけますと幸いです。

1. ニュートン法とは

関数 f(x) が与えられた時、その導関数 f'(x) を用いて、f(x) = 0 となる解 x を反復計算によって求める方法を、ニュートン法 といいます。
ある適当な値 x_0 から計算をはじめ、次の漸化式を繰り返し計算することによって、真の解 α の近似値を求めます。

x_{n + 1} = x_n - \frac{f(x_n)}{f'(x_n)}

ここで、$n = 0, 1, 2, ... $ とする。
この反復計算を実行し、決められた条件が満たされた時に計算を停止し、その時の値 x_n を真の解 α の近似解とする。

具体例として、\sqrt{5} の値を求めてみましょう。

\sqrt{5} の値を求める

\sqrt{5} は、

x^2 - 5 = 0

の解の一つであるので、関数 f(x)

f(x) = x^2 - 5 ...(1)

と定義し、方程式 f(x) = 0 の解をニュートン法で求める。

初期値を

x_0 = 5 ...(2)

とする。f(x) の微分は

f'(x) = 2x ...(3)

であるので、第1ステップ後の値は、(1)(2)(3)から

x_1 = x_0 - \frac{f(x_0)}{f'(x_0)} \\ = 5 - \frac{5^2 - 5}{2 * 5} \\ = 3 ...(4)

となる。第2ステップの解は、(1)(3)(4)より、

x_2 = x_1 - \frac{f(x_1)}{f'(x_1)} = 3 - \frac{3^2 - 5}{2 * 3} = 2.333333333 ...(5)

となる。第3ステップの値は、(1)(3)(5)より、

x_3 = x_2 - \frac{f(x_2)}{f'(x_2)} = 2.333333333 - \frac{(2.333333333)^2 - 5}{2 * 2.333333333} = 2.238095238095238 ...(6)

となる。第4ステップの値は、(1)(3)(6)より、

x_4 = x_3 - \frac{f(x_3)}{f'(x_3)} = 2.238095238095238 - \frac{(2.238095238095238)^2 - 5}{2 * 2.238095238095238} = 2.2360688956433634 ...(7)

となる。第5ステップの値は、(1)(3)(7)より、

x_5 = x_4 - \frac{f(x_4)}{f'(x_4)} = 2.2360688956433634 - \frac{(2.2360688956433634)^2 - 5}{2 * 2.2360688956433634} = 2.236067977499978 ...(8)

となる。
真の値が

\sqrt{5} = 2.23606797749979

であるので、第5ステップの値は、真の値と13桁ほど一致する。
ここで、反復計算を停止する条件を、各ステップ間の変化量の絶対値が 10^{-5} 以下となった時とすると、

|x_5 - x_4| = |2.236067977499978 - 2.2360688956433634| = 9.18 * 10^{-7}

であるから、x_5x_4 の差の絶対値が 10^{-5} 以下になる。よって、計算を第5ステップで停止とする。

2. ニュートン法を実装してみる

前章で紹介した漸化式

x_{n + 1} = x_n - \frac{f(x_n)}{f'(x_n)}

を変形します。ある数値 k の平方根を求める場合、f(x) = x^2 - k となるので

x_{n + 1} = x_n - \frac{x_n^2 - k}{2*x_k} = \frac{1}{2} * (x_n + \frac{k}{x_n})

となります。この式をベースに、ニュートン法を実装しましょう。

newton.py
x = 5.0

while True:
	x2 = x - (x * x - 2) / (x * 2)

	if abs(x2 - x) < 0.0001:
		break

	x = x2

print(x)

比較的シンプルな実装に抑えることができました。


記事内に誤謬等ございましたら、修正いたします。その際は、ご連絡いただけますと幸いです。

それでは、また。

GitHubで編集を提案

Discussion