MastHead

２分検索法 (Binary Search) で解を求める方法


方程式を解く場合に、単純に解が求まらないときがあります。　 それでも、「変数が大きくなるにつれて、従属変数がいつも大きくなる、またはいつも小さくなる」 ということが分かっていて、途中で大きくなったり小さくなったりのような変化がないことが分かっているときは、パソコンで数行の簡単なプログラムを書いてやると解くことができます。
強力な方法ですので非常に利用価値があります。　おぼえると便利です。

まず抽象的な概念から説明します。　
右のグラフは X の２次式で Y が従属変数です。　何かの事前の知識によって、X が正の領域では Y は常に右上がりであることが分かっているとします。　このとき、Y が 0 になるような X を求めることにします。　（従属変数がある値になるような、元の変数を求める）　この例では２次式ですから、平方根の計算で求まってしまいますが、もっと複雑な式の時も応用ができますから、ちょっと我慢を。
概算によって、答えは 0 と 100 の間にあることが判明しているとします。　min を 0、max を 100 とします。
min と max の真ん中 X=50 を使って Y を計算します。　大きすぎたら、max=X として、元に戻り、min と max の真ん中をとって Y を計算します。
もし小さすぎたら min=X としてループします。

こうやっていくと、 X はどんどんもっともらしい値に近づいていきます。

この計算のたびに Y が 0 からどれだけ乖離しているかが分かりますから、その誤差が十分小さく（たとえば0.00000000001 以下に） なったらループから抜け出すことにします。

これを「２分検索アルゴリズム」といいます。　エクセルのユーザ関数を使って書くと次のようになります。
この関数は、エクセルのセルの中で =BSexample(0) として呼び出せば、X の 0 から 100 の範囲にある Y=0 になるような X の値を返します。　実際にやると、0.724745 となります。
Function BSexample(y As Double) As Double  '関数の定義
xmin = 0#: xmax = 100#                      '使用する値の初期化

Loop1:
   x = (xmin + xmax) / 2                    'min と max の真ん中を取る
   ycalculated = 2# * (x + 0.5#) ^ 2 - 3#   'Y を計算する
   If Abs(ycalculated - y) < 0.0000000000001 Then GoTo a10
                      '与えられた y (=0) からの誤差を見て、小さければ抜け出る
   If ycalculated - y > 0 Then xmax = x Else xmin = x
                      'でなければ min か max を置き換え
GoTo Loop1                                   'ループに戻る

a10:
BSexample = x                                '答えを返す
--------------
もうひとつ単純な例を取り上げて解説します。　前回の記事でダイオードの端子電圧を計算したときに使いました。
下の回路ですが、1.5708V の信号で RL=5.6k を通してダイオードに給電すると、ダイオードの電圧はいくらになるでしょうか。
一般的には 0.6V を使う様に教えられますが、前の記事では 1mV まで正確に出さないと歪み率が下がりません。　ですからきっちり計算したくなります。　抵抗の方はリニアですから何とでもなるのですが、ダイオードの方は指数関数ですから単純には計算できません。



原理は簡単です。
上のグラフを見てください。　要は２つの線の交点を求めるわけです。
それで、任意の電圧の時のダイオードの線の電流値から抵抗の線の電流を引いた、差分を見てみます。　

横軸が 1V の時は非常に大きい正の値です。　0V の時は負の値です。　ですから、その２つの間にあることが分かります。
ならば、２つの真ん中を取って、0.5V を見ます。　負の値です。　ならば、1V と 0.5V の間（0.75V）を調べます。　 正の値になります。　ならば、0.75V と 0.5V の真ん中を取って、、、、、、

ということを繰り返していきます。　そのとき、差の量（誤差）を監視していて、誤差がとても小さくなったらループから抜け出すことにします。　 監視する対象は縦軸の誤差量でなく、横軸で、前回からの変化量や上下の差分を使っても良いでしょう。

また、プログラムの間違いがあったときは無限にループしてハングすることがありますから、ループ回数を数えていて、とても大きく （たとえば 10000 回に） なったら、明らかにエラーと分かる値を返してやることにします。

具体例を出しましょう。
エクセルのマクロには２つの記述方法があります。　一つはいわゆるマクロです。　今ひとつはユーザ関数といわれる方法です。　 どちらでもいいですが、セルに組み込んでしまうときはユーザ関数の方が使いやすいでしょう。　 ただし、複数の仮引数は使えませんので、たくさんの答えを持ち帰りたいなら SUB のマクロを使います。

=VBE2SA942(1.5708) という式をセルの中に書いてやると、そのセルは 0.581326 という値になります。
関数のプログラムは次の通りです。
Function VBE2SA942(E As Double) As Double

C = 38.48321418: I0 = 0.000000000000034: RL = 5600#
vmax = 1: vmin = 0: Count = 0

Loop1:
   v = (vmax + vmin) / 2
   Id = I0 * (Exp(v * C) - 1)
   Ir = (E - v) / RL
   If Abs(Id - Ir) < 0.000000000001 Then GoTo a10
   If Id > Ir Then vmax = v Else vmin = v
   Count = Count + 1
   If Count > 10000 Then VBE2SA942 = 9999: GoTo a20
GoTo Loop1

a10:
VBE2SA942 = v

a20:
End Function

それを使っている様子です。　どの場所に書けばいいかが分かるようにしておきました。



ついでですが、上の回路はコンプレッサになっていますので、その入出力特性を出しておきます。　下の図の青い線です。　ひとつ前の記事では、飽和しかかりの曲線部分をうまく使ってサインカーブにフィットさせています。



［余談］
２分検索法は、２つに分けたどちらかに解が存在する、ということを使っているわけですが、必ずしも 「２で」 割らねばならないというわけではありません。　1.5 でも 3 でも 10 でもいいわけです。
私が 34 才のころだと思いますが、それに気づいて、この分割がいくらが効率的なのかを調べたことがあります。　まだ PC のクロックが 1MHz 位だった頃の話です。　結論的には 3 の辺りにあることが分かりました。　そのときは 3.5 が良かったのです。　米国と日本の数学博士の方数人に相談して問題を投げかけてみますと、どなたからも e （ネピア数：自然対数の底）だという理論的な答えが帰ってきました。
ですから、非常に複雑で、計算に時間のかかる関数なら、3 にするといいでしょう。　しかしプログラムが長くなって間違いやすくなりますから、実用的には 2 を使うべきだと思います。　速度差はせいぜい 1.5 倍ぐらいしかありません。
20 歳の頃はまだ２分検索法というのを知らなくて（というか、プログラミングアルゴリズムという学問分野がまだ成熟していなかった）、自分で勝手に 10 分割して 10 進数で１桁ずつ下げていくやり方をしていました。　その方法だと、どこまでが有効数字かということが明白です。 10 進数というものが人間の作ったわがままな数体系だということを理解するよりずっと前のことでした。

［余談２］
ダイオードの電圧を求めるときに、初期値を 1V と 0V から始めました。　理論的には、本来は 1.24V と 0V から始めます。　これは、シリコンのバンドギャップが 1.24V だからです。　 しかし、実用上は 1V から始めて問題ありません。

2015/1/30 追記：
計算機のスピードが当時に比べて数千倍の速さになっています。　そこで、ネピア数が最適である、ということを検証しようと実験をしたところ、ちょっと違う結果が出ました。

0～1 区間の乱数を発生させて、その値を見つけようとしています。　10000 回の平均を出しています。
この実験では 2 が最適になりました。　さて困った。　2 か e か、どちらが本当でしょう。

この問題をゆっくり考えてみました。　私は今は 2 の方が正解だと思います。
まず、1.5 で割る場合(C1.5） と 3 で割る場合(C3)を比較します。
C1.5 では、1.5 の逆数の 0.667 を境にして、その下に答えがあるかを調べます。　無ければ上 1/3 にあることが分かります。　1/3 区間にあるかどうかを調べていることと同じです。
C3 では 3 の逆数 0.333 を境にして、それ以下（1/3 区間）に答えがあるかを調べます。　　なので、 C1.5 と C3 は同等の行為です。
3 の逆数 0.333 と 1.5 の逆数 0.667 とは、0.5 を挟んで対称です。
同様に、2 の逆数 (0.5） を対称軸として、左右にある対の値の逆数による検索はすべて同等になります。

続いて、それでは、その対称図形が上に凸か、下に凸かを考察します。
1 で割ったときは永遠に答えは見つかりません。　また無限大で割ると、無限小区間を調べようとして答えはありません。　従って、下に凸であることが明白です。　　
また、1 を超え無限大未満の区間では特異点はありませんからなめらかなカーブです。

上(1)から下りて来る量と、下(0)から上がっていく量が同じなら同等であるから、真ん中 1/2 で落ち合う、ということです。
（正しいかなあ？）


［余談３］
上のプログラムでは GoTo というインストラクションを使っています。　 本当は、構造化したプログラミング手法ではあまり推奨されないやり方です。
   Do While delta > 0.000000001
     [式]
   Loop
の方がスマートです。