Numerical Recipes in Cの三重対角行列の方程式を解く関数をExcel VBAに移植

ちょっと必要になって三重対角行列の計算をするtridagをVBAに移植してみた。

行列は

A=[ b1 c1 0 .......                0 ]

   [a2   b2 c2 0 ......             0]

   [0    a3 b3 c3 0 ...            0]

   [ 0 ....       aN-1  bN-1  cN-1]

   [0 ..................         aN     bN  ]

のような形としてAx=bを計算する。こんな感じの関数。

Sub tridag(a() As Double, b() As Double, c() As Double, r() As Double, u() As Double, n As Long)
    Dim j As Long
    Dim bet As Double
    Dim gam() As Double
    ReDim gam(n)
   
    If (b(1) = 0#) Then
        MsgBox "Error 1 in tridag"
    End If
   
    bet = b(1)
    u(1) = r(1) / bet
    For j = 2 To n
        gam(j) = c(j - 1) / bet
        bet = b(j) - a(j) * gam(j)
        If (bet = 0#) Then
            MsgBox "Error 2 in tridag"
        End If
        u(j) = (r(j) - a(j) * u(j - 1)) / bet
    Next j
    For j = n - 1 To 1 Step -1
        u(j) = u(j) - gam(j + 1) * u(j + 1)
    Next j
   

End Sub

さて、テストするメインルーチンをこんな感じで書いてみた。

Option Explicit

Private Sub CommandButton1_Click()
    Dim n As Long
    Dim a() As Double, b() As Double, c() As Double
    Dim r() As Double, u() As Double
   
    n = 5
    ReDim a(n), b(n), c(n), r(n), u(n)
   
    b(1) = -5#
    b(2) = -2#
    b(3) = -2#
    b(4) = -1#
    b(5) = -8#
   
    a(2) = 3#
    a(3) = 7#
    a(4) = -2#
    a(5) = 5#
   
    c(1) = 2#
    c(2) = -1#
    c(3) = -4#
    c(4) = 9#
   
    r(1) = -1#
    r(2) = -4#
    r(3) = 0#
    r(4) = 1#
    r(5) = 2#
   
    Call tridag(a, b, c, r, u, n)
   
    MsgBox u(1) & "," & u(2) & "," & u(3) & "," & u(4) & "," & u(5)

End Sub

答えは、1,2,3,2,1でOK。この例は、”C#で学ぶ偏微分方程式の数値解法”に出てたもの。

普通のLU分解で計算する例はこちら。このときはRedimという便利なものを知らなかったので、書きなおす必要があるなあ。

https://sci.tea-nifty.com/blog/2010/07/excel-vbanumeri.html

さて、なぜこの関数を移植したかというと（続く）。