求根アルゴリズム (4) - wetchのブログ

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

少しずつ修正しながらまとめていこう。

まずは標準モジュールに反復計算する本体の関数を書く。

Function IterativeCalculation(f As IMathFunc, ite As IIterationStrategy) As Double
    Const MAX_STEP As Long = 10000
    Const EPSILON As Double = 0.000001

    Dim x() As Double
    x = f.InitialValue
    
    Dim i As Long
    For i = 0 To MAX_STEP
      x = ite.OneCalculation(x, f)
      If Abs(f.Substitute(x(1))) < EPSILON Then Exit For
    Next i
    
    IterativeCalculation = x(1)
End Function

次にクラスモジュールに、数学関数を表すインターフェースクラスIMathFuncを作る。

'定数パラメータのプロパティ。
'複数あれば配列だが1つだけならただのDoubleなので型指定はしないでおく。
Property Get Parameters() 
End Property

'引数のxの型指定はしないでおく。上と同じ理由。
Property Let Parameters(x)
End Property

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

Function InitialValue() As Double()
End Function

同じくクラスモジュールに、反復計算法を表すインターフェースクラスIIterationStrategyを作る。

Function OneCalculation(x() As Double, f As IMathFunc) As Double()
End Function

反復計算法の実装を2つ作る。1つ目は二分法。クラス名をBinaryMethodとする。

Implements IIterationStrategy

Function IIterationStrategy_OneCalculation(x() As Double, f As IMathFunc) As Double()
    Dim result(1) As Double
    result(1) = (x(0) + x(1)) / 2
    
    If f.Substitute(x(0)) * f.Substitute(result(1)) > 0 Then
      result(0) = x(1)
    End If
    
    IIterationStrategy_OneCalculation = result
End Function

同様に2つ目の反復計算法として割線法を実装。名前はSecantMethod。

Implements IIterationStrategy

Function IIterationStrategy_OneCalculation(x() As Double, f As IMathFunc) As Double()
    Dim fx(1) As Double
    fx(0) = f.Substitute(x(0))
    fx(1) = f.Substitute(x(1))
    
    Dim result(1) As Double
    result(0) = x(1)
    result(1) = (x(0) * fx(1) - x(1) * fx(0)) / (fx(1) - fx(0))
    
    IIterationStrategy_OneCalculation = result
End Function

ここまでが準備。以上のモジュールは具体的な問題が変わっても書き変えない。

つまりここ以降は具体的な問題に応じて書き変えていくモジュールになる。
まずは数学関数IMathFuncの実装。名前はたとえばConcreteFunctionとする。

Implements IMathFunc

Private b_ As Double
Private c_ As Double

Property Get IMathFunc_Parameters()
    IMathFunc_Parameters = Array(b_, c_)
End Property

Property Let IMathFunc_Parameters(para)
    b_ = para(0)
    c_ = para(1)
End Property

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

Function IMathFunc_InitialValue() As Double()
    IMathFunc_InitialValue = Array(-2, 1) 'ここに初期値を書く。
End Function

最後に標準モジュールで呼び出しもととなるmainプロシージャを作る。

Sub main()
    Dim f As IMathFunc
    Set f = New ConcreteFunctionOfX '具体的な関数が変わるときはここを書き変える。
    
    '定数パラメータをここで定義する。
    'ConcreteFunctionクラスで定義した配列の順序と合うように注意。
    f.parameters = Array(5, 7)
    
    Dim ite As IIterationStrategy
    Set ite = New BinaryMethod '必要に応じてSecantMethodに書き変える。
   
    Dim x As Double
    x = IterativeCalculation(f:=f, ite:=ite)
    
    Debug.Print "x = " & x
End Sub

とりあえず当初の目的通りには動くものができたかと。
でも改修点はちらほらあるなあ。

main関数の準備の部分が煩雑。
標準モジュールで定数パラメータparametersを設定してるのは、エクセルのワークシートから値を取り込んでいるから。でも初期値の設定はクラスモジュールなので記述場所が離れすぎか？　全部クラスモジュールにしたほうがいいか？
反復解法をもっと実装していこう。以下が参考になる。

求根アルゴリズム二分法 / 割線法 [1/3] - サブロウ丸
 求根アルゴリズム挟み撃ち法 / イリノイ法 / アンダーソン・ビョーク法 [2/3] - サブロウ丸
 求根アルゴリズムデッカー法ブレント法 [3/3] - サブロウ丸
 【Python】 newton法(非線形計画) - サブロウ丸

そもそもの問題は $f(x)=0$ となる $x$ を求めるものだけど、右辺に非ゼロの目標値を与えて $f(x)=a$ となる $x$ を求めるようにしたほうが使い勝手がいいかもしれない。
収束判定条件をもっと合理的なものにしたい。IterativeCalculationでは $|f(x_1)| < 10 ^ {-6}$ としているけど、左辺は一般に有次元な物理量なので無次元量と比較するのは非合理的。
IMathFuncの中に初期値に関する情報を入れるのは、クラスに機能を詰め込みすぎか？
エラー処理を入れていく必要がある。
IterativeCalculationメソッドを標準モジュールに書くと邪魔だからこれもクラスモジュールに移したい。