ニュートン法をExcel VBAで実現してみた

ニュートン法またはニュートン・ラフソン法は、数値計算的に方程式の解を求める手段です。Excelでマクロを組んでみました。

ニュートン法の原理
ニュートン法を使う上での注意事項
1. 増分について
2. 初期値について
Excelマクロ
まとめ

ニュートン法の原理

y = f(x)のように、任意の入力xに対し、関数値yが定まることがニュートン法を使える条件です。このとき、y = 0になるxの値を求めます。関数は微分可能であれば、ブラックボックスで構いません。

もし、y = 5の解を求めたければ、w = y – 5 という関数wを考えて、w = 0の解をニュートン法で求めることに帰着できます。

ニュートン法

さて、例として、y = f(x) = x²– 1という関数を取り上げ、y = 0の解を求めてみましょう。

初期値を(x, y) = (2.5, 5.25)としましょう。この点に接線を引きます。

ここで、元の関数の値が分かっていれば、微係数はy’=2xです。
でも、ブラックボックスだとした場合や微係数が求められない場合は、微係数の定義を用いて近似計算します。

増分を仮に0.01とすると、

$y'\cong\frac{f(2.51)-f(2.5)}{0.01}=\frac{5.3001-5.25}{0.01}=5.01$

と計算されます。確かに、2×2.5に近い値ですね。次に、この微係数を用いて接線を引くと、

となり、この接線とx軸の交点を求めるとy = 0として、

$\begin{eqnarray*} 0-5.25=5.01(x-2.5) \\ x=-\frac{5.25}{5.01}+2.5\cong1.452... \end{eqnarray*}$

となります。

初期値の2.5に比べて、x^2-1 = 0の解の一つであるx = 1に近づいていますね。この手順を何度か繰り返せば、解に収束するという訳です。

ここで、マクロを組むために、式を一般化しておきましょう。

まず、微係数の近似値は、

$y'_0\cong\frac{f\left(x_0+h\right)-f\left(x_0\right)}{x_0+h-x_0}=\frac{f\left(x_1\right)-f\left(x_0\right)}{x_1-x_0}$

です。ただし、初期の状態を(x, f(x)) = (x₀, y₀)とし、(x + h, f(x + h)) = (x₁, y₁)と置きました。
接線は、

$y-y_0=\frac{f\left(x_1\right)-f\left(x_0\right)}{x_1-x_0}\left(x-x_0\right)$

となるので、x軸との交点のx座標は、y = 0として、

$x=x_0-y_0\frac{x_1-x_0}{f\left(x_1\right)-f\left(x_0\right)}$

となります。

ニュートン法を使う上での注意事項

増分について

微係数を求めるわけですから、理論的には増分は無限小であるべきです。

しかし、Excelで無限小を実現することはできません。また、有効数字や丸め誤差の影響も考えないといけません。よって、事例ごとに手入力するか、何かしらのアルゴリズムを使って増分を決める必要があります。

後述するマクロは、初期値の100万分の1を増分としています。

初期値について

y = x² – 1 のグラフ

お気づきの方もいらっしゃると思いますが、x² – 1 = 0の解は±1であり、２つ存在します。初期値をx > 0にとれば、x = 1に収束しますし、x < 0にとればx = -1に収束します。また、初期値がx = 0だと、接線がx軸と平行になるため、収束しません。
すなわち、初期値によって収束値が異なるので、初期値の選び方が大切です。

余談ですが、僕は、アナログ回路技術者です。最適解を得るためには、どうしてもシミュレーションが必要です。最近のシミュレータは、よくできていて、適当な回路を仮定し、ゴールを決めると、それらしい解を出してくれます。よって、少し経験を積むと、シミュレーションだけに頼って、回路を設計することができます。でも、本当のあるべき姿は、手作業で回路構成とパラメータの概算値を決めておいて、シミュレータは最後の微調整だけに使うことだと考えており、これこそが、高性能の回路設計のコツだと信じています。

Excelマクロ

お待たせしました。Excelマクロを紹介します。

まず、新しいシートを開き、「ニュートン法.xlsm」などとして保存します。ファイルの種類を「Excelマクロ有効ブック(*.xlsm)」にしないとマクロを保存することができませんのでご注意ください。

A1、A2セルにそれぞれ「x=」「y=」と入力します。
次にB1、B2セルにそれぞれ「x」と「y」の名前を付けます。
B1セルを選択し、左上のB1と表示されている▼印の付いた小窓に「x」とタイプし、Enterを押します。これで、B1セルがxという文字で参照できるようになります。次に、B2セルを選択し、同様に「y」という名前を付けます。マクロからこれらのセルにアクセスする際には名前を使うので、これをしていないと、マクロが動きません。

さらに、B1、B2セルにそれぞれ「2.5」「=x^2-1」と入力します。前者が初期値、後者が関数です。B2は、入力が終わると「5.25」と表示されるはずです。

次に、「開発」タブから「Visual Basic」をクリックします。
「開発」タブが見当たらない場合や、マクロの使い方が初めての方は、こちらの記事を参照してください。

左のプロジェクトエクスプローラで、「VBAProject（ニュートン法.xlsm）」を右クリックし、「挿入」「標準モジュール」を選びます。そこに開いたエディタに下記のコードをコピー&ペーストしてください。

'-----------------------------------------------------------------------
' Summary: ニュートン法を実行する
'-----------------------------------------------------------------------
Public Sub Newton()
    ' ─── 計算の過程を表示する列の準備 ───
    ' 計算の過程を表示する列を定義する
    Const cstColumnX As String = "D"    ' Xを表示する列
    Const cstColumnY As String = "E"    ' Yを表示する列
    ' 計算の過程を表示する列をクリアする
    Columns(cstColumnX & ":" & cstColumnY).Clear
    ' 1行目にXとYのタイトルを表示する
    Const cstRowTitle As Integer = 1
    Range(cstColumnX & cstRowTitle) = "x"
    Range(cstColumnY & cstRowTitle) = "y"
    
    ' ─── ニュートン法を実行する ───
    Dim i       As Integer  ' ループ変数
    Dim dblH    As Double   ' Xの増分
    ' 回数の上限を決めてループ
    Const cstRowMax As Integer = 100
    For i = cstRowTitle + 1 To cstRowMax
        Dim dblX0 As Double ' 初期のX値
        Dim dblY0 As Double ' 初期のY値
        ' 初期のX値を保持する
        dblX0 = Range("x").Value
        ' 初期のY値を保持する
        dblY0 = Range("y").Value
        ' 計算の過程を表示する列にもコピーしておく
        Range(cstColumnX & i).Value = Range("x").Value
        Range(cstColumnY & i).Value = Range("y").Value
        ' Xの増分を決める(初期値の100万分の1)
        If dblH = 0 Then
            dblH = dblX0 / 1000000#
            If dblH = 0 Then
                dblH = 0.000001
            End If
        End If
        ' 初期値からXが増分だけ増えた地点の座標を求める
        Dim dblX1 As Double
        Dim dblY1 As Double
        ' Xを増分だけ増加させる
        Range("x").Value = dblX0 + dblH
        ' 座標を保持する
        dblX1 = Range("x").Value
        dblY1 = Range("y").Value
        ' 微係数が0だと値を求めることができないので終了する
        If dblY1 - dblY0 = 0 Then
            Call MsgBox("傾きが0なので他の初期値を入れてください", _
                    vbOKOnly + vbCritical, "エラー")
            Exit Sub
        End If
        ' オーバフローに備えてエラーの処理方法をセット
        On Error Resume Next
        ' 接戦とx軸の交点を求める
        Range("x").Value = dblX0 - dblY0 * (dblX1 - dblX0) / (dblY1 - dblY0)
        ' エラーが発生していたら内容を表示して終了する
        If Err.Number <> 0 Then
            Call MsgBox(Err.Description, vbOKOnly + vbCritical, "エラー")
            Exit Sub
        End If
        ' エラー処理を解除する
        On Error GoTo 0
        ' 値に変化がなければ終了したとみなす
        If Range("x").Value = dblX0 And Range("y").Value = dblY0 Then
            Exit For
        End If
    Next i
End Sub