こんにちは、はっこんです。
今回は数値計算法の一つである、二分法をpythonで作ってみたいと思います。
原理
二分法とは方程式の解を出力する方法で、数値計算法の中で最も単純で理解しやすいものだと思います。
原理を説明すると、
- 下限aと上限bをきめる
- aとbの中間点cを求める
- f(c)の符号がf(a)と同じならばcをaに置き換え、f(c)の符号がf(b)と同じならばcをbに置き換える
- 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位までは正確な値が出力されるということになりますね。
ここまで理解すれば、かなり単純であることが分かると思います。
このような方程式を求める数値計算方は、二分法の他にもニュートン法などがあげられますね。
まぁ今日のところはここまでにしときましょう。
コメント