Math＆Julia #08｜二分法とニュートン・ラフソン法 » 平方根

はじめに

方程式を数値的に解くための代表的な手法として「二分法（bisection）」 と 「ニュートン・ラフソン法（Newton-Raphson）」があります。本稿ではこの 2つの手法を使って平方根を求める方程式を解きます。

平方根を求める方程式

ある数 $a>0$ の平方根を求めるには「 $2$乗すると $a$ になる数 $x$ 」を探します。この条件を関数 $f(x)$ で表すと次のようになります。

$$\tag{1.1}f(x)=x^2-a$$

知りたいのは $f(x)=0$（ $x$ をいくつにすれば $f(x)$ の値が $0$ になるのか？）なので、次の方程式を解くことになります。

$$\tag{1.2}{f(x)=x^2−a=0}$$

実例

$a=2$ の平方根を探した結果、解 $x=\sqrt{2}$ を見つけたとすると次のように方程式が満たされます。

$$\tag{1.3}f\big(\sqrt{2}\,\big)=\big(\sqrt{2}\,\big)^2-2=0$$

中間値の定理

関数 $f(x)$ が区間 $[\,x_{\text{low}},\,x_{\text{high}}\,]$ で連続であり、かつ次の式が真であるとき、この区間内には $f(x)=0$ となる解$x$が必ず存在します。

$$\tag{2.1}f(x_{\text{low}}) \cdot f(x_{\text{high}}) < 0$$

実例

$a=2$ の平方根を探しているとします。式$(1.1)$と式$(2.1)$を用いて、区間 $[\,1,\,2\,]$ に解が存在するか調べます。

$$\tag{2.2}f(1)\cdot f(2)=(1^2-2)\cdot(2^2-2)=-1\cdot 2=-2<0$$

この結果より、区間 $[\,1,\,2\,]$ に解が存在します（ $x=\sqrt{2}$ ）。

二分法（bisection）

二分法の基本的なアイデアは、中間値の定理に基づいて解が存在する区間を半分ずつ狭めていくというものです。

図 1　二分法で解（√2）を求める

初期区間の設定

初期区間 $[\,x_{\text{low}},\,x_{\text{high}}\,]$ は次のように設定します。これは中間値の定理を満たしており適切です。

$$\tag{3.1}[\,x_{\text{low}},\,x_{\text{high}}\,]=[\,0,\,\text{max}(1,\,a)\,]$$

中点の計算

区間の中点を求める方法は直感的に分かる通りです。この中点 $x_{\text{mid}}$ が最終的に近似解になります。

$$\tag{3.2}x_{\text{mid}}=\dfrac{x_{\text{low}}+x_{\text{high}}}{2}$$

区間の更新

$f(x_{\text{mid}})$ の符号を調べると区間の更新方法が分かります。次の式を用いて区間を狭めていくことによって、中点 $x_{\text{mid}}$ が解に近づいて行きます。

$$\tag{3.3} \begin{cases} \enspace x_{\text{low}}=x_{\text{mid}}&\text{if}&f(x_{\text{mid}})<0\\[4pt] \enspace x_{\text{high}}=x_{\text{mid}}&\text{if}&f(x_{\text{mid}})>0 \end{cases}$$

収束判定

次の収束判定が真になったら $x_{\text{mid}}$ を近似解として反復計算を打ち切ります。

$$\tag{3.4}|\,x_{\text{high}}-x_{\text{low}}\,|<ε$$

$$または$$

$$\tag{3.5}|\,f(x_{\text{mid}})\,|<ε$$

• $ε$は許容誤差

実例

図 1では式$(3.2)$を用いて中点を計算しています。また、図には明示していませんが区間の更新には式$(3.3)$を用いています。

ニュートン・ラフソン法（Newton–Raphson）

ニュートン・ラフソン法の基本的なアイデアは、関数の接線と $x$軸の交点を利用して解に近づいて行くというものです。

図 2　ニュートン・ラフソン法で解（√2）を求める

ニュートン・ラフソン法の漸化式を導出

点$(x_n,\,f(x_n))$ における接線の方程式は、直線の方程式 $y-f(x_n)=m(x-x_n)$ より次のようになります。

$$\tag{4.1}y-f(x_n)=f'(x_n)(x-x_n)$$

この接線が $x$軸と交わる点は、$x=x_{n+1},\:y=0$ となる点です。

$$\tag{4.2}0-f(x_n)=f'(x_n)(x_{n+1}-x_n)$$

これを $x_{n+1}$ について解くとニュートン・ラフソン法の漸化式が導出されます。

$$\begin{align} \tag{4.3}&-\dfrac{f(x_n)}{f'(x_n)}=x_{n+1}-x_n\\[8pt] \tag{4.4}&x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)} \end{align}$$

平方根の漸化式を導出

今回はニュートン・ラフソン法を使って平方根を求めたいのでした。まずは、平方根の関数とその導関数を用意します。

$$\begin{align} \tag{5.1}&f(x)=x^2−a\\[4pt] \tag{5.2}&f'(x)=2x \end{align}$$

式$(5.1)$と式$(5.2)$をニュートン・ラフソン法の漸化式$(4.4)$に代入します。

$$\tag{5.3} x_{n+1}=x_n-\dfrac{x_n^2-a}{2x_n}=\dfrac{2x_n^2-(x_n^2-a)}{2x_n}=\dfrac{x_n^2+a}{2x_n}$$

これをさらに変形すると平方根の漸化式が導出されます。

$$\tag{5.4}x_{n+1}=\dfrac{1}{2}\,\bigg(x_n+\dfrac{a}{x_n}\bigg)$$

これは結果的にはバビロニア法（古代バビロニアにおける平方根の計算方法）として知られている式です。

平方根の漸化式に初期値を設定

平方根の漸化式の初期値 $x_1$ は次のように設定します。

$$\tag{5.5} x_1=\begin{cases} \enspace a/2&\text{if}&a>1\\ \enspace 1&\text{if}&a≤1 \end{cases}$$

図 2では $x_1=a$ に従って初期化していますが、式$(5.5)$の $x_1=a/2$ に従うと、厳密解から極端に離れていない初期値を設定できます。

収束判定

次の収束判定が真になったら $x_{n+1}$ を近似解として反復計算を打ち切ります。

$$\tag{5.6}|\,x_{n+1}-x_n\,|<ε$$

• $ε$は許容誤差

実例

実例として図 2の $x_1,\,x_2,\,x_3$ を式$(5.4)$を用いて手計算すると次のようになります。

$$\begin{aligned} &x_1=a=2\\[4pt] &x_2=0.5\,(2+2/2)=1.5\\[4pt] &x_3=0.5\,(1.5+2/1.5)=1.4166\cdots \end{aligned}$$

実装

二分法とニュートン・ラフソン法の 2つの手法について、平方根の数値計算を実装します。両者とも反復計算を行いますが、その収束速度は大きく異なるので、これを視覚的に確認できるようにします。また、計算結果（近似解）については、高精度解と比較することにより精度を把握できるようにします。

Julia

using CairoMakie
using Printf

function sqrt_bisect(a, ε; maxiter=60)  # 平方根（二分法）
    f(x)   = x ^ 2 - a                  # 平方根の関数（区間の設定で使用）
    lo, hi = 0.0, max(1.0, a)           # 区間の初期値
    xs     = Float64[]                  # 途中経過（中点）の記録用
    for i in 1:maxiter
        x = (lo + hi) / 2               # 区間の中点を求める
        push!(xs, x)                    # 中点を記録
        hi - lo < ε && break            # 許容誤差内であれば早期打ち切りへ
        f(x) > 0 ? hi = x : lo = x      # 区間を狭めて再試行へ
    end
    xs
end

function sqrt_newton(a, ε; maxiter=20)  # 平方根（ニュートン・ラフソン法）
    f(x) = 0.5 * (x + a / x)            # 平方根の漸化式
    x    = a > 1 ? a/2 : 1.0            # 漸化式の初期値
    xs   = Float64[x]                   # 途中経過（反復値）の記録用
    for i in 1:maxiter
        x_new = f(x)                    # 平方根の反復計算
        push!(xs, x_new)                # 反復値を記録
        abs(x_new - x) < ε && break     # 許容誤差内であれば早期打ち切りへ
        x = x_new                       # 次の反復計算へ
    end
    xs
end

function raw_convergence(vals, val_true, a, title)  # 平方根の近似解
    val = vals[end]                     # 近似解を取得
    err = abs(val - val_true)           # 誤差を計算
    acc = (1 - err/val_true) * 100      # 精度を計算
    println("$title\n", "-"^40)
    @printf "被開平数　　　　 : %s\n"        a
    @printf "反復回数　　　　 : %d\n"        length(vals)
    @printf "平方根‐近似解　 : %.17f\n"     val
    @printf "平方根‐高精度解 : %.77f\n"     val_true
    @printf "誤差　　　　　　 : %.17f\n"     err
    @printf "精度　　　　　　 : %.17f%%\n\n" acc
end

function vis_convergence(vals_bisect, vals_newton, val_true)  # 収束速度の比較
    errs_bisect = abs.(vals_bisect .- val_true)  # 誤差（二分法）
    errs_newton = abs.(vals_newton .- val_true)  # 誤差（ニュートン・ラフソン法）
    fig = Figure()
    ax  = Axis(fig[1,1], yscale=log10, xlabel="Iteration", ylabel="Error")
    scatterlines!(ax, errs_bisect, color=:blue, marker=:circle,  label="Bisection")
    scatterlines!(ax, errs_newton, color=:red,  marker=:diamond, label="Newton-Raphson")
    axislegend(ax)
    display(fig)
end    

ハンズオン

Julia

a = 2                                   # 被開平数
ε = 1e-16                               # 許容誤差
val_true    = sqrt(big(a))              # 平方根（高精度解）
vals_bisect = sqrt_bisect(a, ε)         # 平方根（二分法）
vals_newton = sqrt_newton(a, ε)         # 平方根（ニュートン・ラフソン法）
vis_convergence(vals_bisect, vals_newton, val_true)  # 収束速度の比較
raw_convergence(vals_bisect, val_true, a, "二分法")
raw_convergence(vals_newton, val_true, a, "ニュートン・ラフソン法") 

図 3　収束分析（√2）

実行結果

二分法
----------------------------------------
被開平数　　　　 : 2
反復回数　　　　 : 60
平方根‐近似解　 : 1.41421356237309492
平方根‐高精度解 : 1.41421356237309504880168872420969807856967187537694807317667973799073247846210
誤差　　　　　　 : 0.00000000000000013
精度　　　　　　 : 99.99999999999999113%

ニュートン・ラフソン法
----------------------------------------
被開平数　　　　 : 2
反復回数　　　　 : 7
平方根‐近似解　 : 1.41421356237309492
平方根‐高精度解 : 1.41421356237309504880168872420969807856967187537694807317667973799073247846210
誤差　　　　　　 : 0.00000000000000013
精度　　　　　　 : 99.99999999999999113%

ニュートン・ラフソン法は二分法と比べて少ない反復回数で急速に収束しているのが分かります。また、両者の近似解を比較すると同等の精度を確保していることが分かります（表示範囲の限りでは同じです）。

本稿の内容だけだとニュートン・ラフソン法の方が一方的に優れているように見えますが、こちらは条件によっては収束しないことがあります。一方の二分法は中間値の定理を満たせば必ず収束するという性質があります。それぞれの性質についての詳細は省きますが、実用面では両者の長所を組み合わせたハイブリッド戦略も有効とされています。