AI&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}$ ）。

二分法

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

図 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)$を用いています。

ニュートン法

ニュートン法の基本的なアイデアは、関数の接線と $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, 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)                    # 中点を記録
        abs(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 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")
    axislegend(ax)
    display(fig)
    fig
end    

# 平方根の近似解と高精度解をテキストで表示
function show_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

実行例

初期設定

2の平方根を求めたいので、被開平数には 2を設定します。

Julia

a = 2                                   # 被開平数
ε = 1e-16                               # 許容誤差
val_true = sqrt(big(a))                 # 平方根（高精度解）

二分法とニュートン法の収束分析

二分法とニュートン法を用いて近似解を計算しますが、先に収束の様子を観察します。

Julia

vals_bisect = sqrt_bisect(a, ε)         # 平方根（二分法）
vals_newton = sqrt_newton(a, ε)         # 平方根（ニュートン法）
fig = vis_convergence(vals_bisect, vals_newton, val_true)  # 収束の比較
save("08-conver.png", fig)

図 3　収束分析（√2）

ニュートン法は二分法と比べて少ない反復回数で急速に収束しているのが分かります。

二分法とニュートン法の近似解

二分法とニュートン法の近似解を比較します。また、近似解と高精度解も比較してみます。

Julia

show_convergence(vals_bisect, val_true, a, "二分法")
show_convergence(vals_newton, val_true, a, "ニュートン法")

実行結果

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

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

二分法とニュートン法の近似解を比較すると同等の精度を確保していることが分かります（表示範囲の限りでは同じです）。近似解と高精度解を比較しても十分な精度を確保していることが分かります。

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

おわりに

2の平方根を求めるという最小限の題材を通じて、二分法とニュートン法について理解することができました。特にニュートン法については勾配降下法に通じる考え方が含まれている点が有益でした。