pythonで遊んでみる〜二分法

知識

こんにちは、はっこんです。

今回は数値計算法の一つである、二分法をpythonで作ってみたいと思います。

原理

二分法とは方程式の解を出力する方法で、数値計算法の中で最も単純で理解しやすいものだと思います。

原理を説明すると、

  1. 下限aと上限bをきめる
  2. aとbの中間点cを求める
  3. f(c)の符号がf(a)と同じならばcをaに置き換え、f(c)の符号がf(b)と同じならばcをbに置き換える
  4. 2,3を繰り返し、f(x)=0の解の値に近づけていく

って感じですね。

探索範囲を半分に狭めていき、解を出力するというわけです。

めっちゃ単純ですよね。

しかし、二分法の弱点は、

解の範囲がある程度分かっていること、

探索範囲に解が1つしかないことが確実、というときにしか使えないことですね。

ということで実装してみましょう!

コード

Scipyを用いて

Scipyという科学技術計算ライブラリのoptimizeというサブモジュールをインポートすることでとても簡単に二分法を利用することができます。

今回は例として、

f(x)=x^2-2

という方程式の解を出力してみます。

コードはこんな感じ。

from scipy import optimize

def f(x):
    return x**2 - 2

print(optimize.bisect(f, 0, 2))

出力結果がこちら、

1.4142135623715149

今回の方程式の解はもちろん

\sqrt{2}

ですが、一般的に覚えられている範囲までは求まっていますね。

(ひとよひとよにひとみごろ)

Scipyを用いればたった数行で二分法を実装することができます。

Scipyを使わずに

もちろんScipyは非常に便利なライブラリですが、Scipyを使わずに二分法を実装することもできます。

プログラミングの知識が浅いボクが言うのもおこがましいですが、

できる限り関数がどのようなコードで構成されているかを知ることも大切なことです。

ということでScipyを使わずに、二分法を実装してみましょう。

こんなかんじ。

def f(x):
    return x**2 - 2

def bisect(f,a,b):
    eps = 0.00001
    print("n\tc\tf(c)")
    n = 1
    while True:
        c = (a + b)/2
        print("{}\t{:.5f}\t{:.5f}".format(n, c, f(c)))
        if f(a)*f(c) < 0:
            b = c
        else:
            a = c
        n += 1
        if abs(f(c)) < eps or abs(a - b) < eps:
            break
    
    print("x = %f" % c)

a = 0
b = 5
bisect(f,a,b)

出力結果がこちら、

n       c       f(c)
1       2.50000 4.25000
2       1.25000 -0.43750
3       1.87500 1.51562
4       1.56250 0.44141
5       1.40625 -0.02246
6       1.48438 0.20337
7       1.44531 0.08893
8       1.42578 0.03285
9       1.41602 0.00510
10      1.41113 -0.00870
11      1.41357 -0.00181
12      1.41479 0.00164
13      1.41418 -0.00008
14      1.41449 0.00078
15      1.41434 0.00035
16      1.41426 0.00013
17      1.41422 0.00003
18      1.41420 -0.00003
19      1.41421 -0.00000
x = 1.414213

だんだんと解に近づいていることが分かります。

ループ処理は、

f(c)<\varepsilon\,(\varepsilon=1\times10^{-5})

もしくは、

|a-b|<\varepsilon\,(\varepsilon=1\times10^{-5})

を満たすまで続きます。

ということは、少数第5位までは正確な値が出力されるということになりますね。

ここまで理解すれば、かなり単純であることが分かると思います。

このような方程式を求める数値計算方は、二分法の他にもニュートン法などがあげられますね。

まぁ今日のところはここまでにしときましょう。

コメント

タイトルとURLをコピーしました