0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

ゆるゆる数値計算Advent Calendar 2018

Day 8

二分法の数理

Last updated at Posted at 2018-12-07

はじめに

この記事では二分法を扱います。数値計算では良く知られた手法で、簡単に実装できるため教科書でも最初の方に出てくるものでしょう。その解説をします。

適用対象

二分法の適用対象は、非線形実数値連続関数です。次の問題を考えます。

区間$I$上で定義された連続関数$f:I \to \mathbb{R}$に関して、$f(c) = 0, c \in I$となる$a$を計算することを考えます。

もちろん一般にこの式は成り立ちません。例えば、 $f(x) = x^2 + 1$としたら実数区間の中でどのように$I$をとっても$f(c) = 0$となることはないでしょう。
逆に$I$の取り方次第では、たくさんあることもあります。例えば、 $f(x) = \cos x$とかなら$f(c) = 0$となる$c$はたくさんありますよね。

中間値の定理

二分法の数学的基礎は中間値の定理です。杉浦1の第1章 実数と連続 §8 中間値の定理によれば

定理8.1(中間値の定理) $\mathbb{R}$の有界閉区間 $I = [a, b]$で連続な実数値関数$f$は、 $f(a)$と$f(b)$の間の任意の実数$\gamma$を値にとる:$f(c) = \gamma$となる$c \in I$が存在する。
定理8.1系1  $f: I = [a, b] \to \mathbb{R}$が連続で、$f(a) < 0,\ f(b) > 0$(または、$f(a) > 0,\ f(b) < 0$)ならば、$f(c) = 0$となる$c \in I$が存在する。

系1はまさに解きたい問題ですよね。つまり「区間が定義されていて、その端の符号が異なっている」というのが二分法の適用条件です。実は二分法はこの証明そのものといった感じです。次の節で証明します。

中間値の定理の証明

定理8.1で $\tilde{f}(x) = f(x) - \gamma$とすれば、系1の仮定を満たすから、系1を証明すれば良い。
ここで$f(a) < 0 < f(b)$としても一般性を失わない。 ここで系列$a_n, b_n, n\in\mathbb{N}$を$a_0 = a, b_0 = b,c_n = \frac{a_n + b_n}{2}$とした時に、

a_{n+1} = a_n, b_{n+1} = c_n\ {\rm if}\ f(c_n) > 0\\
a_{n+1} = c_n, b_{n+1} = b_n\ {\rm otherwise}

と定義し、区間$I_n$を$I_n = [a_n, b_n]$で定義する。この時、区間の系列$I_n$は有界閉区間の単調減少列で、その長さは$\frac{b - a}{2^n}$なので$0$に収束する。故に、区間縮小法より$I_n \to \{c\}, a_n \to c, b_n \to c$となる。一方、全ての$n$に対して、

f(a_n) \leq 0 \leq f(b_n)

が成り立つ。$n=0$の時には仮定から、$f(a_0) < 0 < f(b_0)$であり、系列の構成法から、

f(a_{n+1}) = f(a_n) \leq 0 \leq f(b_{n+1}) = f(c_n)\ {\rm if}\ f(c_n) > 0\\
f(a_{n+1}) = f(c_n) \leq 0 \leq f(b_{n+1}) = f(b_n)\ {\rm otherwise}

が成り立つため、帰納法から証明できる。ここで$f$は連続であったから、$f(c) \leq 0 \leq f(c)$より$f(c) = 0$である。

解いてみる

さて、せっかくなのでコードを書いて実際に問題を解いてみましょう。今回のコードはpythonで書きました。

def bs_method(f, start, end, N=100):
    if f(start) * f(end) >= 0:
        raise Exception('not satisfying the solvable condition.')
    # 初期化
    a, b = start, end
    s_fa = f(a) > 0
    # 以下をN回(初期値100)繰り返す
    for _ in range(N):
        # 中間の値の取得
        med = (a + b) / 2
        f_med = f(med)
        # もし中間の値が解ならそれを返す
        if f_med == 0:
            return med
        else:
            s_fmed = f_med > 0
        # 符号が同じ方を交換する
        if s_fmed == s_fa:
            a = med
        else:
            b = med
    return med
bs_method(lambda x: 2 *  x * x  - 1, 0, 5) ==> 0.7071067811865475

このように二分法では繰り返しの回数を指定する必要があります。
この回数は精度に依存します。また、最初に与えた区間の大きさも制度に依存します。
$k$回の繰り返しをしたとして、精度はおおよそ$(b - a)/2^k$になります。
必要な精度で使うと良いです。

おまけ

さらに微分が計算できる場合、ニュートン法という方法があります。
こちらは収束が基本的には良いですが、解の周りが平らに近い場合、数値が発散しやすいという特徴を持っています。
個人的には、2階微分があり、convexかconcaveかがわかっているならニュートン法やその改良法、そうでないなら二分法という使い分けをします。

  1. 杉浦光夫 著, 「解析入門I」 東京大学出版会, 1980

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?