二分法を用いて関数の根を求める

1変数連続関数の根(値が0になる点)を求める、二分法というアルゴリズムとそれをScipyで実装する方法を紹介します。

アルゴリズムの内容

二分法というのは中間値の定理をベースとした求根アルゴリズムです。アイデアは非常に単純で、連続関数$f$に対して、$f(x_1)$と$f(x_2)$の符号が異なるように、$x_1, x_2$を選びます。この時点で、中間値の定理より区間$(x_1, x_2)$に根があることがわかりますのでさらに細かくみていきます。次は、$x_1, x_2$の中点$x_M = \frac{x_1+x_2}{2}$を取り、$f(x_M)$の符号を調べます。$f(x_M)$と$f(x_1)$の符号が同じであれば、$x_1$を$x_M$で置き換え、逆に$f(x_M)$と$f(x_2)$の符号が同じであれば、$x_2$を$x_M$で置き換えると、区間の幅が半分になった$(x_1, x_2)$が得られますが根はこの中にあることがわかります。これを繰り返すと、根が存在する範囲を狭めていくことができ、$f(x_M)$の絶対値が0になるか、もしくは十分0に近づいたらその値を数値的に求めた根とします。

以上が、一般的な二分法のアルゴリズムの説明です。ただし、後に紹介するSciPyではどうやら区間が十分狭くなったかイレーション回数が上限に達したか等の基準でループを打ち切り、$f(x)$の値を確認していないようです。

二分法のメリットとデメリット

方法が単純でわかりやすい、というのが個人的に感じている１番のメリットです。

また、連続関数であれば使えるため、ニュートン法などのアルゴリズムのように元の関数の微分を必要とせず、微分が難しい関数や微分不可能な関数でも使えます。

また、根が存在しうる区間を狭めながら探索するため、最初の区間の幅と繰り返し回数により、結果の精度を保証できることも大きな利点です。

逆にデメリットとしては、ニュートン法等と比較して収束が遅いこととか、初期値として関数が異符号になる2点を探して与える必要があること、もしその区間に複数の根が存在した場合にどれに収束するか不確定なことなどが挙げられます。

ただし僕の経験上では、ある程度根の目処がついていたり、単調な関数に対して使う場面が多くこれらのデメリットを深刻に感じることは少ないです。

SciPyによる実装

SciPyではscipy.optimizeというモジュールで実装されています。bisectという専用メソッドを使うか、root_scalarという汎用的なメソッドで(method=’bisect’)を指定して使うことになります。

参考:
– scipy.optimize.bisect — SciPy v1.11.4 Manual
– root_scalar(method=’bisect’) — SciPy v1.11.4 Manual

試しに、$\sin$関数の 3と4の間にある根を探させて見せましょう。既知の通りそれは円周率$\pi$になるはずです。

from scipy import optimize
import numpy as np


root1 = optimize.bisect(np.sin, 3, 4)
print(root1)
# 3.1415926535901235

root_result = optimize.root_scalar(np.sin, bracket=[3, 4], method='bisect')
# 結果は複数の情報を含むRootResults形式で戻る。
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 3.1415926535901235
"""
print(type(root_result))
# class 'scipy.optimize._zeros_py.RootResults'

# 根の値へのアクセス方法
print(root_result.root)
# 3.1415926535901235

いかにも円周率ぽい結果が得られましたね。root_scalarの方では収束したことを示すフラグや、イテレーション回数なども得られています。

失敗事例1. 区間の両端での関数の値が同符号の場合

二分法は初期設定を誤ってると失敗するのでその場合のSciPyの挙動も見ておきましょう。失敗パターンの一つは、最初に指定した区間の両端で符号が一致していた場合です。もちろん関数の形によってはその区間内に根がある可能性もあるのですが、存在は保証されなくなります。

また$\sin$関数で、その間に根が存在しない区間$(1, 2)$と、実は両端で同符号だけど根が存在する区間$(1, 7)$でやってみましょう。bisectとroot_scalarで全く同じエラーメッセージ出るのでbisectの方だけ載せます。

try:
    optimize.bisect(np.sin, 1, 2)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

try:
    optimize.bisect(np.sin, 1, 7)
except Exception as e:
    print(e)
# f(a) and f(b) must have different signs

はい、$f(a)$と$f(b)$は違う符号にせよとのことでした。根が存在しない1個目の例はさておき、2個目の例は根は区間内に2個存在するのですが探さずにエラーになりました。

失敗事例2. 連続関数ではなかった場合

もう一つ失敗するのは、関数が連続関数ではないケースです。

例えば$\tan$の、$\frac{\pi}{2}$近辺の挙動で見てみましょう。数学的に厳密な話をすると、$\frac{\pi}{2}$では$\tan$は定義されないので、$\tan$は数学的には連続関数(定義域内のすべての点で連続)なのですが、数値計算的には不連続と考えた方が都合が良いです。

話が脇に逸れたので実例の話に移ります。実は、bisectメソッドは結果を返して来ちゃうんですよね。そしてそれが全然根ではないということも見ておきましょう。

root = optimize.bisect(np.tan, 1, 2)
# pi/2に近い結果が得られている
print(root)
# 1.5707963267941523

# 元の関数に代入した結果は全く0に近くない
np.tan(root)
# 1343445450736.3804

root_scalarの方だったら、結果のフラグ等もあるのでアラート等あげてくれるのかと期待したのですがそういう機能はなさそうです。

root_result = optimize.root_scalar(np.tan, bracket=[1, 2], method='bisect')
print(root_result)
"""
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 1.5707963267941523
"""
np.tan(root_result.root)
# 1343445450736.3804

要するに、SciPyに渡す関数が本当に連続関数であるかどうかは利用者が責任を持たないといけないということです。また、結果が本当に根なのかどうかは代入して確認した方が良いでしょう。

まとめ

以上が二分法の説明とSciPyで利用する方法、その注意点でした。