ニュートン法 その2

前回お話したニュートン法を，EXCELで利用するためのVBAモジュールについてお話します
EXCELにはゴールシーク等数値解を探索するツールもあるんですがちょっと使いにくいので・・・（'A`)
newton.xls
【書式】
= newton(expr, a)
expr : 数式（文字列）
a : 初期値（Double）

• 非線形方程式 expr = 0 の解をニュートン法で一次探索し，得られた数値解を返します
• exprはxの関数として表現して下さい
• 数値解は初期値aに依存します
• 所定の探索回数(imax=1000)で収束しなかった場合は 10^99 を返します
• 数式を直接記入する場合は，数式を""で囲んで下さい（他セルから参照させる場合は不要です）

【入力例】

F列 : 解きたい非線形方程式を入力（ = 0 は省略）
a列 : 初期値を入力
x列 : "=newton(F列, a列)"を入力，返値は数値解
F(x)列 : "=f(F列, x列), "を入力，返値は数値解における残差
入力例の3つとも残差がほぼ0となっているため，十分な精度の数値解が得られていることが解ります

※標準モジュールに下記を追加し，マクロを有効にして下さい

Function newton(expr As String, a As Double) As Double

    Dim i, imax As Integer
    Dim eps, Xnew, X As Double

    imax = 1000
    eps = 10 ^ -6

    X = a

    For i = 1 To imax

        Xnew = X - f(expr, X) / df(expr, X)
        If (Abs(Xnew - X) < eps * a) Then Exit For
        X = Xnew

    Next i
    
    If (i = imax) Then Xnew = 10 ^ 99

    newton = Xnew

End Function


Function f(expr, X As Double) As Double

    f = Evaluate(Replace(expr, "x", X))

End Function


Function df(expr, X As Double) As Double

    dx = X * 10 ^ -4
    df = (f(expr, X + dx) - f(expr, X)) / dx

End Function