ミュラー法とは？わかりやすく解説

ミュラー法（ミュラーほう、ミューラー法、英語: Muller's method）は、求根アルゴリズムであり、f(x) = 0 の形式の方程式を解く数値解析の手法である。 1956年にDavid E. Muller（英語版）によって初めて発表された。

ミュラー法は、セカント法（割線法）の2次漸化式に類似した3次漸化式に従って進行する。セカント法では、関数 f のグラフ上の、直前2回の近似値に対応する2点を通る直線を構築し、その直線の根を新たな近似値として用いることを反復する。一方、ミュラー法では、直前3回の近似値に対応する3点を使用する。これら3点を通る放物線を構築し、その放物線の根を新たな近似値として用いることを反復する。

導出

ミュラー法では、3つの根の初期近似値 $x_{0},x_{1},x_{2}$ を用いる。そして $(x_{0},f(x_{0}))$ , $(x_{1},f(x_{1}))$ , $(x_{2},f(x_{2}))$ を通る放物線とx軸の交点を考慮することによって、新たな近似値 $x_{3}$ を決定する。

$(x_{0},f(x_{0}))$ , $(x_{1},f(x_{1}))$ , $(x_{2},f(x_{2}))$ を通る2次多項式

P(x)=a(x-x_{2})^{2}+b(x-x_{2})+c

(1)

について、差分

$h_{0}=x_{1}-x_{0},\quad h_{1}=x_{2}-x_{1}$

および $\delta _{0}={\frac {f(x_{1})-f(x_{0})}{h_{0}}},\quad \delta _{1}={\frac {f(x_{2})-f(x_{1})}{h_{1}}}$

を定義する。 3つの点 $(x_{0},f(x_{0}))$ , $(x_{1},f(x_{1}))$ , $(x_{2},f(x_{2}))$ をそれぞれ式(1)に代入し、 $a,b,c$ について解くと、

$a={\frac {\delta _{1}-\delta _{0}}{h_{1}+h_{0}}},\quad b=ah_{1}+\delta _{1},\quad c=f(x_{2})$

を得る。次に、この2次方程式を(1)に適用して $x_{3}$ を求めると、

x_{3}-x_{2}={\frac {-2c}{b\pm {\sqrt {b^{2}-4ac}}}}

となる。次の反復で $x_{2}$ に最も近くなるように、根号の手前の符号は $b$ の符号に一致させる。すると、

x_{3}=x_{2}-{\frac {2c}{b+\operatorname {sign} (b){\sqrt {b^{2}-4ac}}}}

となる。 $x_{3}$ を決定したら、この手順を繰り返す。なお、分母に含まれる平方根の式のため、前回までの反復ではすべて実数であっても、次の反復では複素数になる可能性がある。これは、セカント法、Sidiの一般化セカント法（英語版）、ニュートン法などのように、実数から開始した場合に反復も実数のままになる他の根探索アルゴリズムとは対照的である。反復が複素数であることは、問題によって利点（複素根を探している場合）または欠点（すべての根が実数であることがわかっている場合）になる可能性がある。

この手法は衝突検出のように、最も近い根ではなく、より小さい根が必要な用途にも、 $\operatorname {sign} (b)$ を $\operatorname {sign} (-a)$ に置き換えることで容易に適用できる：

x_{3}=x_{2}-{\frac {2c}{b-\operatorname {sign} (a){\sqrt {b^{2}-4ac}}}}.

収束速度

良好な関数の場合、ミュラー法の収束次数（英語版）は約1.839であり、トリボナッチ定数と一致する。これは、黄金比と一致するセカント法の約1.618や、ニュートン法の2と比較することができる。つまり反復1回あたりの進捗が、セカント法はミュラー法よりも小さく、ニュートン法の方が大きいということである。

より正確には、ξ が f の単根を表す場合（つまり f(ξ) = 0 かつ f'(ξ) ≠ 0）、f は3回連続微分可能であり、初期推定値 x₀, x₁, x₂ が ξ に十分近いとすると、反復は

\lim _{k\to \infty }{\frac {|x_{k}-\xi |}{|x_{k-1}-\xi |^{\mu }}}=\left|{\frac {f'''(\xi )}{6f'(\xi )}}\right|^{(\mu -1)/2}

を満たす。ただし μ ≈ 1.84 は、トリボナッチ定数の定義式 $x^{3}-x^{2}-x-1=0$ の正の解である。

一般化と関連手法

ミュラー法では、各反復において、最後に得られた3点 f(x_k-1), f(x_k-2), f(x_k-3) に放物線、つまり2次多項式を当てはめる。これを一般化して、k 回目の反復で、最後の m+1 点に次数 m の多項式 p_k,m(x) をあてはめることができる。この記法では、放物線y_k はp_k,2と表される。次数 m は1以上でなければならない。新たな近似値 x_k は、p_k,m の根の1つ、つまり p_k,m(x)=0 の解の1つになる。 m=1 とするとセカント法になり、m=2 とすると Muller 法になる。

ミュラーは、このように生成された数列 {x_k} が、次数 μ_m で根 ξ に収束することを計算した。ただし、μ_m は $x^{m+1}-x^{m}-x^{m-1}-\dots -x-1=0$ の正の解である。

m が無限大に近づくにつれて、この方程式の正の解は2に近づく。しかし m>2 の場合、この方法は m=1 または m=2 の場合よりも格段に難しくなる。これは、3次以上の多項式の根を求めることが難しいためである。さらに m>2 の場合、p_k,m のどの根を新たな近似値 x_k として選ぶべきかについて、明確な基準がないように見える、との問題もある。

これらの困難は、Sidiの一般化セカント法（英語版）によって克服される。この方法でも多項式 p_k,m を用いるが、p_k,m(x)=0 を解く代わりに、x_k-1 における p_k,m の導関数を用いて新たな近似値 x_k を計算する。

実装例

以下は、ミュラー法をPythonプログラミング言語で実装したものである。このプログラムは、根の初期推定値3つ、期待精度（小数点以下の桁数）、最大反復回数をパラメータとして受け取る。このプログラムは、関数 $f (x) = x 2 - 612$ の根を求める目的に適用される。

Pythonでの実装例

from cmath import sqrt  # Use the complex sqrt as we may generate complex numbers

def muller(x0: int, x1: int, x2: int, decimal_places: int, maximum_iterations: int):
    func = lambda x: (x ** 2) - 612
    iteration_counter = 0
    iterates = [x0, x1, x2]
    solution_found = False

    while not solution_found and iteration_counter < maximum_iterations:
        final_index = len(iterates) - 1

        h0 = iterates[final_index - 1] - iterates[final_index - 2]
        h1 = iterates[final_index] - iterates[final_index - 1]

        f_x0 = func(iterates[final_index - 2])
        f_x1 = func(iterates[final_index - 1])
        f_x2 = func(iterates[final_index])

        delta0 = (f_x1 - f_x0) / h0
        delta1 = (f_x2 - f_x1) / h1

        coeff_a = (delta1 - delta0) / (h1 + h0)
        coeff_b = coeff_a * h1 + delta1
        coeff_c = f_x2

        sqrt_delta = sqrt(pow(coeff_b, 2) - 4 * coeff_a * coeff_c)
        denominators = [coeff_b - sqrt_delta, coeff_b + sqrt_delta]

        # Take the higher-magnitude denominator
        next_iterate = iterates[final_index] - (2 * coeff_c) / max(denominators, key=abs)

        iterates.append(next_iterate)
        solution_found = abs(func(next_iterate)) < pow(10, -decimal_places)
        iteration_counter = iteration_counter + 1

    if solution_found:
        print("Solution found: {}".format(next_iterate))
    else:
        print("No solution found.")

muller(10, 20, 30, 9, 20)  # Solution found: (24.73863375370596+0j)

脚注

参考文献

Atkinson, Kendall E. (1989). An Introduction to Numerical Analysis, 2nd edition, Section 2.4. John Wiley & Sons, New York. ISBN 0-471-50023-2.
Burden, R. L. and Faires, J. D. Numerical Analysis, 4th edition, pages 77ff.
Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). “Section 9.5.2. Muller's Method”. Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8

ミュラー法とは？わかりやすく解説