数値的逆ラプラス変換(vba-code)

Rem 数値的逆ラプラス変換filt06(t) Hosonoのアルゴリズム----------------------
Rem atstp:分母 btstp:分子---ラプラス関数値の算出に必要なパラメータをセル参照
Rem 作例 F(s) = ht(8次多項式)/gt(8次多項式) : 有理関数
Public Function filt06(time As Double, atstp As Variant, btstp As Variant, nts As Integer) As Double
Dim ss(500) As Complex
Dim x1L(500) As Complex
Dim res(500) As Double
Dim aa As Double
Dim n As Integer
Dim cnt As Double

Rem パラメータ関連・複素数変数
Dim resigma As Double
Dim ared(50, 2) As Double
Dim bred(50, 2) As Double
Dim pi As Double
Dim a01(50) As Complex
Dim b01(50) As Complex
Dim ht(50) As Complex
Dim gt(50) As Complex
Dim ht01 As Complex
Dim ht02 As Complex
Dim ht03 As Complex
Dim ht04 As Complex
Dim gt01 As Complex
Dim gt02 As Complex
Dim gt03 As Complex
Dim gt04 As Complex

Rem 各種定数
pi = 3.14159265358979

Rem パラメータの取り込み・複素数変換
Rem 指定範囲から順番にデータ用セル値を取り出す(次数の高い順a0-a8,b0-b8)
For Each myCell In atstp
      ared(cnt, 0) = myCell.Value
      a01(cnt) = ToComplex(ared(cnt, 0), 0)
      cnt = cnt + 1
Next
cnt = 0
For Each myCell In btstp
      bred(cnt, 0) = myCell.Value
      b01(cnt) = ToComplex(bred(cnt, 0), 0)
      cnt = cnt + 1
Next

Rem 高速ラプラス変換(細野法)
If time > 0 Then
aa = 8

For n = 1 To nts + 10
         ss(n) = ToComplex(aa / time, (n - 0.5) / time * pi)

        Rem ラプラス変換関数 Cadd Csub Cmul Cdiv ToComplex F(ss(n))算出 
        gt(7) = ss(n)
        gt(6) = Cmul(gt(7), ss(n))
        gt(5) = Cmul(gt(6), ss(n))
        gt(4) = Cmul(gt(5), ss(n))
        gt(3) = Cmul(gt(4), ss(n))
        gt(2) = Cmul(gt(3), ss(n))
        gt(1) = Cmul(gt(2), ss(n))
        gt(0) = Cmul(gt(1), ss(n))

        gt01 = a01(8)
        gt01 = Cadd(gt01, Cmul(gt(7), a01(7)))
        gt01 = Cadd(gt01, Cmul(gt(6), a01(6)))
        gt01 = Cadd(gt01, Cmul(gt(5), a01(5)))
        gt01 = Cadd(gt01, Cmul(gt(4), a01(4)))
        gt01 = Cadd(gt01, Cmul(gt(3), a01(3)))
        gt01 = Cadd(gt01, Cmul(gt(2), a01(2)))
        gt01 = Cadd(gt01, Cmul(gt(1), a01(1)))
        gt01 = Cadd(gt01, Cmul(gt(0), a01(0)))

        ht(7) = ss(n)
        ht(6) = Cmul(ht(7), ss(n))
        ht(5) = Cmul(ht(6), ss(n))
        ht(4) = Cmul(ht(5), ss(n))
        ht(3) = Cmul(ht(4), ss(n))
        ht(2) = Cmul(ht(3), ss(n))
        ht(1) = Cmul(ht(2), ss(n))
        ht(0) = Cmul(ht(1), ss(n))

        ht01 = b01(8)
        ht01 = Cadd(ht01, Cmul(ht(7), b01(7)))
        ht01 = Cadd(ht01, Cmul(ht(6), b01(6)))
        ht01 = Cadd(ht01, Cmul(ht(5), b01(5)))
        ht01 = Cadd(ht01, Cmul(ht(4), b01(4)))
        ht01 = Cadd(ht01, Cmul(ht(3), b01(3)))
        ht01 = Cadd(ht01, Cmul(ht(2), b01(2)))
        ht01 = Cadd(ht01, Cmul(ht(1), b01(1)))
        ht01 = Cadd(ht01, Cmul(ht(0), b01(0)))

        x1L(n) = Cdiv(ht01, gt01)
        Rem ラプラス変換関数 Cadd Csub Cmul Cdiv ToComplex

        res(n) = Imz(x1L(n)) * (-1) ^ n
Next

Rem nts項分加算
resigma = 0
For n = 1 To nts
        resigma = resigma + res(n)
Next

Rem オイラー変換-10項分追加
resigma = resigma + (res(nts + 1) * 1023 + res(nts + 2) * 1013 + res(nts + 3) * 968 + res(nts + 4) * 848 + res(nts + 5) * 638 ) / 2 ^ 10
resigma = resigma + (res(nts + 6) * 386 + res(nts + 7) * 176 + res(nts + 8) * 56 + res(nts + 9) * 11 + res(nts + 10)) / 2 ^ 10

Rem 変換値出力
filt06 = resigma * Exp(aa) / time
Else
Rem t=0の関数値:初期値の出力
filt06 = 0
End If
End Function

 有理関数の数値的逆ラプラス変換vbaのコードです。ユーザー定義関数の形にしてあります。解説は後程です。