目的

二分法により一変数間数 f(x)=0 の解を求める。
R には uniroot という関数もあり，これはおもしろい（有効な）使い方ができる。

使用法

bisection(func, lower, upper, ndiv=50, epsilon=1e-14, max.iteration=500)

引数

fun           関数定義
lower         解を探索する範囲の下限値
upper         解を探索する範囲の上限値
ndiv          解を探索する範囲の区分数（50）
epsilon       許容誤差（1e-14）
max.rotation　収束計算上限界数（500）

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/bisection.R", encoding="euc-jp")

# 分法により一変数間数 f(x)=0 の解を求める
bisection <- function(       func,                                   # 関数　たとえば function(x) (x+6.7)*(x-3.4) のようなもの
                        lower, upper,                           # 解を求める区間
                        ndiv=50,                                # 区間を細区分する個数
                        epsilon=1e-14,                          # 解の許容誤差
                        max.iteration=500)                      # 反復回数の上限
{
        printf <- function(fmt, ...) cat(sprintf(fmt, ...))  # 書式付きの print
        bisec2 <- function(func, lower, upper)                       # 区間の再区分について解を探索する下請け関数
        {
                yl <- func(lower)                            # 区間の左端における関数値
                yu <- func(upper)                            # 区間の右端における関数値
                if (yl*yu > 0) {                                # 同じ符号のとき，この区間には解がない（区間の設定不良）
                        return(1)                               # 戻り値として 1 を返す（利用しなくてもよい）
                }
                for (i in 1:max.iteration) {                    # 繰り返し上限まで収束計算を続ける
                        mid <- (lower+upper)/2                       # そのときの区間の中点の値
                        ym <- func(mid)                              # 中点における関数値
                        if (abs(ym) < epsilon) {             # たまたま解になったら結果を書き出す
                                printf("ans=%g\n", mid)
                                return(0)
                        }
                        else if (yu*ym > 0) {                   # 区間の右端での関数値と符号が同じなら，
                                upper <- mid                 # 中点を区間の右端にする
                                yu <- ym                     # 新しい区間の右端での関数値を設定する
                        }
                        else {                                  # 区間の左端での関数値と符号が同じなら，
                                lower <- mid                 # 中点を区間の左端にする
                                yl <- ym                     # 新しい区間の左端での関数値を設定する
                        }
                        if (abs(upper-lower) < epsilon) {    # 区間幅が誤差範囲内になったら中点を解とする
                                printf("ans=%g\n", (lower+upper)/2)
                                return(0)
                        }
                }
        }

        x <- seq(lower, upper, length.out=ndiv)                      # 区間を細区分して，
        for (i in 1:(ndiv-1)) {                                 # それぞれの細区分に対して，
                bisec2(func, x[i], x[i+1])                      # 二分法により解を求める下請け関数を呼ぶ
        }
}


使用例

bisection(function(x) (x+6.7)*(x-3.4)*(x-sqrt(2)), -10, 10, ndiv=300)		# 3 次方程式の解

出力結果例

> bisection(function(x) (x+6.7)*(x-3.4)*(x-sqrt(2)), -10, 10, ndiv=300)
ans=-6.7
ans=1.41421
ans=3.4