多変数関数の求根アルゴリズムを切り替える

3つの実引数をとる関数 $f(x_1, x_2, x_3)$ について、次のそれぞれには独立している3つの問題を考える：

3つのうち $x_1$ 以外には既知の定数 $X_2, X_3$ を与えて、 $x_1$ を動かすことで方程式 $f(x_1, X_2, X_3)=0$ を数値的に解きたい。ただし $0\leq x_1\leq 1$
$x_2$ 以外には既知の定数 $X_1, X_3$ を与えて、 $x_2$ を動かすことで方程式 $f(X_1, x_2, X_3)=0$ を解きたい。ただし $0\leq x_2\leq 10$
$x_3$ 以外には既知の定数 $X_1, X_2$ を与えて、 $x_3$ を動かすことで方程式 $f(X_1, X_2, x_3)=0$ を解きたい。ただし $-2\leq x_3\leq 2$

ただし $f$ の単調性は仮定し、解も存在することは保証されているとする。

求根アルゴリズムとして二分法を採用することとする。
二分法 - Wikipedia

すると問題1についてはこんな感じになる。

※Dimを省略しているとか数値をハードコーディングしてるとかコード中に日本語を使ってるとか、いろいろ適当に書いている点は目をつぶってください。

Function x1を求める(X2, X3)
  x1_left = 0 '左側の初期値
  x1_right = 1 '右側の初期値

  For i = 0 To 10000
    x1_mid = (x1_left + x1_right) / 2
    
    If f(x1_mid, X2, X3) * f(x1_left, X2, X3) > 0 Then
      x1_left = x1_mid
    Else
      x1_right = x1_mid
    End If

    '収束判定
    If Abs(f(x1_mid, X2, X3)) < 10 ^ -6 Then Exit For
  Next i

  x1を求める = x1_mid
End Function

問題2，3についてもほぼ同じようなコードになる。あえて書いてみます。

Function x2を求める(X1, X3)
  x2_left = 0
  x2_right = 10

  For i = 0 To 10000
    x2_mid = (x2_left + x2_right) / 2
    
    If f(X1, x2_mid, X3) * f(X1, x2_left, X3) > 0 Then
      x2_left = x2_mid
    Else
      x2_right = x2_mid
    End If

    '収束判定
    If Abs(f(X1, x2_mid, X3)) < 10 ^ -6 Then Exit For
  Next i

  x2を求める = x2_mid
End Function

Function x3を求める(X1, X2)
  x3_left = -2
  x3_right = 2

  For i = 0 To 10000
    x3_mid = (x3_left + x3_right) / 2
    
    If f(X1, X2, x3_mid) * f(X1, X2, x3_left) > 0 Then
      x3_left = x3_mid
    Else
      x3_right = x3_mid
    End If

    '収束判定
    If Abs(f(X1, X2, x3_mid)) < 10 ^ -6 Then Exit For
  Next i

  x3を求める = x3_mid
End Function

こんなの誰が見たって、コピペコードが多すぎるのでまとめて1つのコードに書きたくなる。

どの変数を変化させるかをString型引数で与えることにして、切り替えをswitchで行うようにする。

Function x_iを求める(変化させる変数 As String, X1, X2, X3)
  Select Case 変化させる変数
    Case "x1"
      y_left = 0
      y_right = 1
    Case "x2"
      y_left = 0
      y_right = 10
    Case "x3"
      y_left = -2
      y_right = 2
  End Select

  For i = 0 To 10000
    y_mid = (y_left + y_right) / 2
    
    Select Case 変化させる変数
      Case "x1"
        f_mid = f(y_mid, X2, X3)
        f_left = f(y_left, X2, X3)
      Case "x2"
        f_mid = f(X1, y_mid, X3)
        f_left = f(X1, y_left, X3)
      Case "x3"
        f_mid = f(X1, X2, y_mid)
        f_left = f(X1, X2, y_left)
    End Select

    if (f_mid * f_left > 0) then
      y_left = y_mid
    Else
      y_right = y_mid
    End If

    '収束判定
    Select Case 変化させる変数
      Case "x1"
        f_mid = f(y_mid, X2, X3)
      Case "x2"
        f_mid = f(X1, y_mid, X3)
      Case "x3"
        f_mid = f(X1, X2, y_mid)
    End Select

    If Abs(f_mid) < 10 ^ -6 Then Exit For
  Next i

  x_iを求める = y_mid
End Function

・・・select caseが鬱陶しすぎる。共通化の恩恵が全然得られていない。なんとかならんものか。

というのをもう5年くらい放置考えていて、この間解決のヒントが見つかった気がしたので書いてみる。なんとここまでが今回の話の前段。

要は各 $x_i$ を選んだ時に初期値と $f$ への代入の仕方だけが違うのだから、そこだけ抜き出せばいい。あとの処理は共通。

標準モジュールから説明すると、次のように書き変える。最初の1か所だけselect caseを使うけど、あとは全部なくせる。

Function x_iを求める(変化させる変数 As String, x1, x2, x3)
  Dim var As IVariable
  Select Case 変化させる変数
    Case "x1"
      Set var = New x1
    Case "x2"
      Set var = New x2
    Case "x3"
      Set var = New x3
  End Select

  y_left = var.初期値left
  y_right = var.初期値right

  For i = 0 To 10000
    y_mid = (y_left + y_right) / 2
    
    If var.fに代入(y_left, x1, x2, x3) * var.fに代入(y_mid, x1, x2, x3) > 0 Then
      y_left = y_mid
    Else
      y_right = y_mid
    End If

    '収束判定
    f_mid = var.fに代入(y_mid, x1, x2, x3)
    If Abs(f_mid) < 10 ^ -6 Then Exit For
  Next i

  x_iを求める = y_mid
End Function

その代わりクラスモジュールを4つ作る。まずはインターフェースIVariable

Function 初期値left()
End Function

Function 初期値right()
End Function

Function fに代入(y, x1, x2, x3)
End Function

あとの3つはIVariableを実装するクラスモジュール。1つ目はx1。

Implements IVariable

Function IVariable_初期値left()
  IVariable_初期値left = 0
End Function

Function IVariable_初期値right()
  IVariable_初期値right = 1
End Function

Function IVariable_fに代入(y, x1, x2, x3)
  IVariable_fに代入 = f(y, x2, x3)
End Function

次に同じくIVariableを実装するクラスモジュールx2。

Implements IVariable

Function IVariable_初期値left()
  IVariable_初期値left = 0
End Function

Function IVariable_初期値right()
  IVariable_初期値right = 10
End Function

Function IVariable_fに代入(y, x1, x2, x3)
  IVariable_fに代入 = f(x1, y, x3)
End Function

最後にクラスモジュールx3

Implements IVariable

Function IVariable_初期値left()
  IVariable_初期値left = -2
End Function

Function IVariable_初期値right()
  IVariable_初期値right = 2
End Function

Function IVariable_fに代入(y, x1, x2, x3)
  IVariable_fに代入 = f(x1, x2, y)
End Function

リファクタリングの本で見たいわゆるStrategyパターンというのを参考に考えた書き方なのだけど、これで合ってるのか・・・？
結局クラスモジュールはコピペで作ってるし、たとえば $f$ が4変数関数になったときに修正する箇所が多すぎる気もする。

あと、二分法を別のアルゴリズム、例えば割線法に変えたいときのことを考えて、forループの中もクラス化したい。

続く

wetchのブログ

他人に見られることを想定していない書き散らかし独習ノート．物理学とかVBAとか．

多変数関数の求根アルゴリズムを切り替える