平方根の値をニュートン法で求めてみる!

Warning: count(): Parameter must be an array or an object that implements Countable in /home/c9962019/public_html/mochinoki-labo.com/wp-content/plugins/rich-table-of-content/functions.php on line 490

こんにちは、MSKです。
今回は平方根の値をニュートン法で求めてみたいと思います。

ニュートン法とは

方程式 $f(x)=0$ の解xを近似的に求めるアルゴリズムです。

ニュートン法

$f$ が解を含むある区間 $[a,b]$ で、次を満たすとする。

2回微分可能
単調増加で下( or 上)に凸
$f$ と $f’$ が計算可能

この時、初期値を $x_0$ として、
$x_{n+1} = x_n – \frac{f(x_n)}{f'(x_n)}$
を繰り返すことで、 $f(x)=0$ の $[a,b]$ 上の解を求めることができる。

注意) 上の解は単調性より、唯一つであることがわかる。

略証といか、軽く解説をしてみます。（数学にあまり興味ない方は飛ばしてください。）

上の状況を数式で整理すると、

$f(a) < 0 < f(b)$
$f'(x) > 0$ $(\forall x \in [a,b])$
$f”(x) > 0$ $(\forall x \in [a,b])$

中間値の定理より、
$\exist c \in (a,b) ~~ s.t.~~ f(c) = 0$
以下、簡単のために $x_0=b$ とします。

$f(b)>0$ より、 $c < b = x_0$ 。
今、 $x_i$ を中心としたテイラー展開により、ある $y_i \in [x_i,x_{i+1}]$ が存在して、

$0=f(c) = f(x_i)+f'(x_i)(c-x_i) + \frac{f”(y_i)}{2}(c-x_i)^2$ ・・・①

を満たします。

$f”(y_i) > 0$ なので、 $\frac{f”(y_i)}{2}(c-x_i)^2 > 0$

よって、上の①から $\frac{f”(y_i)}{2}(c-x_i)^2$ を引いたものは0以下。

すなわち、 $0\geq f(x_i)+f'(x_i)(c-x_i)$

式変形を行って、
$f(x_i)+f'(x_i)c – f'(x_i)x_i \leq 0$
$f'(x_i)c \leq f'(x_i)x_i-f(x_i)$
$f'(x_i)$ で両辺を割って、
$c \leq x_i-\frac{f(x_i)}{f'(x_i)} = x_{i+1}$

よって、 $x_i$ は $[c,b]$ の点列。
また、 $f$ は単調増加なので、 $f(x_i) \geq f(c) = 0$ 。

$f'(x_i)>0$ だったので、
$x_i > x_i – \frac{f(x_i)}{f'(x_i)} = x_{i+1}$
であることが分かる。

よって、数列 $\{x_i\}^{\infty}_0$ は単調減少列。
下に有界な単調減少列は収束する。

最後に $f(x_{\infty})=0$ を示せば、この方法により解を求められることが分かる。（つまりは $c = x_{\infty}$ ）

$x_{i+1} = x_i – \frac{f(x_i)}{f'(x_i)}$ において、
$i \to \infty$ とすると、

$\frac{f(x_\infty)}{f'(x_\infty)} = 0$

仮定より、 $f'(x_\infty)>0$ なので、
$f(x_\infty)=0$ が得られた。

ニュートン法で平方根を計算してみる

mathモジュールを使って、計算しておく

まずは計算があっているか確認するためにmathモジュールを使って、平方根を計算してみます。
今回は $\sqrt{19}$ を計算します。

import math

print(math.sqrt(19))

結果は4.358898943540674という数値になります。

Pythonで計算してみる

今、 $x^2-19$ の正の解として $\sqrt{19}$ を計算したいので、範囲を $[0,5]$ として、初期値を5とします。

# 19の平方根を計算するための式を定義
def polynomial(x: float) -> float:
    return x**2 - 19

#上の式の微分
def polynomial_derivative(x: float) -> float:
    return 2*x

# ニュートン法に出てくる式を計算する
def calc_newton_method_sequence(x: float) -> float:
    return x - (polynomial(x) / polynomial_derivative(x))

# ニュートン法を計算
def calc_newton_method(x: float) -> float:
    val = x
    loop_cnt = 0
    while loop_cnt <= 10:
        val = calc_newton_method_sequence(val)
        loop_cnt += 1
    return val

if __name__ == '__main__':
    print("root 19 =", format(calc_newton_method(5.0))

これを実行するとroot 19 = 4.358898943540674と表示されると思います。
上で計算した値と一致しています。

Rustで計算してみる

ついでにRustでも計算してみます。

// 19の平方根を計算するための式を定義
fn polynomial(x: f64) -> f64 {
    x*x - 19.0
}

// 上の式の微分
fn polynomial_derivative(x: f64) -> f64 {
    2.0 * x
}

// ニュートン法に出てくる式を計算する
fn calc_newton_method_sequence(x: f64) -> f64 {
    x - polynomial(x) / polynomial_derivative(x)
}

// ニュートン法を計算
fn calc_newton_method(x: f64) -> f64 {
    let mut val: f64 = x;
    let mut loop_cnt: u8 = 0;
    while loop_cnt <= 10 {
        val = calc_newton_method_sequence(val);
        loop_cnt += 1;
    }
    val
}

fn main() {
    println!("root 19 = {}",calc_newton_method(5.0));
}

Pythonの時と同じようにroot 19 = 4.358898943540674と表示されます。

最後に

ニュートン法を紹介して、平方根を計算してみました。
ニュートン法は $f(x)=0$ の解xを近似的に求めるアルゴリズムでした。
数列の計算を繰り返し行うことで解に限りなく近づくことができます。

実際に計算すると大変ですが、単純な計算はコンピューターの力を借りると簡単に計算できます。その実例として、PythonとRustでニュートン法による計算を行ってみました。

最後までご覧いただき、ありがとうございます。
「平方根の値をニュートン法で求めてみる!」でした。