wetchのブログ

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

求根アルゴリズム (2)

前回の続き。
https://wetch.hatenablog.com/entry/2019/10/02/212051

前回は1つの多変数関数を持ってきて変化させる変数を切り替えられるようにしたが、仕様を変更&追加する。

  1. 任意に1変数関数f(x; a)を与えて切り替えられるようにしたい。ただしaは定数パラメータベクトル。
  2. 二分法の部分を割線法に切り替えられるようにしたい。割線法 - Wikipedia
  3. mainプロシージャからこれらの切り替えを行わせるようにする。反復法を行うFunction部分には変更を入れなくてもいいようにしたい。

まずは1番と2番から対応しよう。このためにIVariableとその実装を修正する。名前も単に数学関数を表しているという意味なのでIMathFuncに変更。

'定数パラメータのプロパティ
Property Get Parameters()
End Property

Property Let Parameters()
End Property

Function Substitute(x) '変化させる変数に特定の値を代入する
End Function

Function 初期値Left()
End Function

Function 初期値Right()
End Function

実装としてはたとえばこんな感じ。変化させたい変数がx1の場合は次のクラスConcreteFunctionOfX1を使う。

Implements IMathFunc
Private parameters_()

Property Get IMathFunc_Parameters() As Double()
  IMathFunc_Parameters = parameters
End Property

Property Let IMathFunc_Parameters(para() As Double)
  parameters_ = para
End Property

Function IMathFunc_Substitute(x)
  b = parameters_(0)
  c = parameters_(1)

  IMathFunc_Substitute = x * (b - c) ^ 2 'ここに具体的な関数形を書く。
End Function

Function IMathFunc_初期値Left()
  IMathFunc_初期値Left = 0 '初期値も具体的な問題が与えられるたびに書き変える。
End Function

Function IMathFunc_初期値Right()
  IMathFunc_初期値Right = 1
End Function

同様にして、変化させたい変数がx2の場合は次のクラスConcreteFunctionOfX2を使う。

Implements IMathFunc
Private parameters_()

Property Get IMathFunc_Parameters() As Double()
  IMathFunc_Parameters = parameters
End Property

Property Let IMathFunc_Parameters(para() As Double)
  parameters_ = para
End Property

Function IMathFunc_Substitute(x)
  a = parameters_(0)
  c = parameters_(1)

  IMathFunc_Substitute = a * (x - c) ^ 2 'ここが変わっていることに注意
End Function

Function IMathFunc_初期値Left()
  IMathFunc_初期値Left = 0
End Function

Function IMathFunc_初期値Right()
  IMathFunc_初期値Right = 10 'ここも変わっていることに注意
End Function

x3を変えたい場合も同様。

次に3番、mainプロシージャからこれらの切り替えを行わせるようにする。

Sub main()
  Dim f As IMathFunc
  Set f = New ConcreteFunctionOfX1 '変化させたい変数に応じて、ConcreteFunctionOfX2とかを使う。

  Dim parameters(0) '定数パラメータを配列parametersにまとめる。
  parameters(0) = X2
  parameters(1) = X3

  x1 = 反復法(f, parameters)
End Sub

Function 反復法(f As IMathFunc, parameters())
  y_left = f.初期値Left
  y_right = f.初期値Right

  For i = 0 To 10000
    y_mid = (y_left + y_right) / 2
    
    If f.Substitute(y_left) * f.Substitute(y_mid) > 0 Then
      y_left = y_mid
    Else
      y_right = y_mid
    End If

    '収束判定
    If Abs(f.Substitute(y_mid)) < 10 ^ -6 Then Exit For
  Next i

  反復法 = y_mid
End Function

具体的な関数形を与えるIMathFuncクラス(とその実装のConcreteFunctionクラス)で初期値を与えるのはなんとなく違う気がしているけど、とりあえず今回はここまで。
今後は初期値を決める部分も、前回は決め打ちで書いていたが実際はもっと複雑な(f(x;a)に応じて変化する)アルゴリズムになるので切り替え可能にしたい。

続く