はとむぎ研究室

大学生の備忘録でした

【数値計算#1】根を求める

Python 数値計算

根ってなんだよ

　ある関数f(x)が0になるときのxの値のこと。
例：f(x)=x^{2}-1の根は\pm{1}である。

二分法

概要

　初期区間[ x_l,x_r]のどちらかを中点\frac{x_1+x_2}{2}に適切に置き換え、どんどん幅を狭めて根を求めていく。
　長所：確実である
　短所：収束は遅い

詳しく説明してみる

　y=x^{2}-2について考えてみる。次のグラフは、y=x^{2}-2 (赤)とy=0 (黒)を図示している。

f:id:heart_mugi:20201001212228j:plain
二分法(1)
　さて、初期値[ x_l,x_r]を[ 1,2]とする。これも図示してみる。
f:id:heart_mugi:20201001212347j:plain
二分法(2)
　次に、x_lとx_rの中点x_m=を求める。中点x_m=\frac{x_l+x_r}{2}=\frac{1+2}{2}=1.5となり、これも図示してみる。
f:id:heart_mugi:20201001212533j:plain
二分法(3)
　ここで、f(x_l)(x_m)とf(x_m)(x_r)について考える。もし、とある2点からなる幅の間に根があるとき、その2点をそれぞれ関数に代入した値どうしの積は負の値になるはず(正\times負=負)である。f(x_l)f(x_m)\lt 0, f(x_m)f(x_r)>0より、[ x_l(=1),x_m(=1.5)]を採用。なので、x_r=x_mと置き換える。
f:id:heart_mugi:20201001213550j:plain
二分法(4)
　以下繰り返し。f(a)f(b)\lt 0かつ|a-b|\lt \epsilonとなったとき、探索を終了する。直前の中点が答えになる。

Pythonソースコード

def bisection(func, xl, xr):
    epsilon = 1.0e-6
    i = 0
    while(abs(xr-xl) > epsilon):
        sl = func(xl)
        sr = func(xr)
        xm = (xl + xr) / 2    # 中点
        sm = func(xm)
        if(sl*sm < 0.0):
            xr = xm
        else:
            xl = xm
        i += 1
        print('iteration:{}, value:{}'.format(i,xm))
    
    return xm

実際に動かしてみる

関数y=x^{2}-2

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


実行

bisection(func, 1, 2)

iteration:1, value:1.5
iteration:2, value:1.25
iteration:3, value:1.375
iteration:4, value:1.4375
iteration:5, value:1.40625
iteration:6, value:1.421875
iteration:7, value:1.4140625
iteration:8, value:1.41796875
iteration:9, value:1.416015625
iteration:10, value:1.4150390625
iteration:11, value:1.41455078125
iteration:12, value:1.414306640625
iteration:13, value:1.4141845703125
iteration:14, value:1.41424560546875
iteration:15, value:1.414215087890625
iteration:16, value:1.4141998291015625
iteration:17, value:1.4142074584960938
iteration:18, value:1.4142112731933594
iteration:19, value:1.4142131805419922
iteration:20, value:1.4142141342163086