こんばんは、皆様、三頌亭です。数値的逆ラプラス変換のGaver-Stehfest 法についての作例をあげておきます。係数はn=14に固定してあります。逆変換の例はもっとも簡単な指数関数をやってみました。細野法と2つを今までの記事で説明してきましたが、算法の解説に終始してきました。実際使用する人の立場として関連の数学の論文読んでも「具体的にどうやるんだよ!!」という人も多いのではないかと思います。計算法・解説マニュアルめいたものが「BASICによる高速ラプラス変換」(絶版)ぐらいしかないのも残念です。決して悪い手法ではないと思うのですが・・・。そんな思いがこれらの記事を書いた理由でもあります。
参照記事
https://kms130.hatenablog.com/entry/2020/12/05/165813
Rem f(t)=exp(-alpha*t) : F(s)=1/(s+alpha)
Public Function Stehfest01(time As Double, alpha As Double) As Double
Dim ind(40) As Double
Dim reg As Double
Dim n As Double
Dim ln As Double
ln = 0.693147180559945
ind(1) = 0.0027777778
ind(2) = -6.4027777778
ind(3) = 924.05
ind(4) = -34597.9277777778
ind(5) = 540321.111111111
ind(6) = -4398346.36666667
ind(7) = 21087591.7777778
ind(8) = -63944913.0444444
ind(9) = 127597579.55
ind(10) = -170137188.083333
ind(11) = 150327467.033333
ind(12) = -84592161.4999999
ind(13) = 27478884.7666666
ind(14) = -3925554.96666666
If time > 0 Then
For n = 1 To 14
Rem ラプラス変換関数記述
reg = reg + 1 / (ln * n / time + alpha) * ind(n)
Rem ラプラス変換関数記述
Next n
Stehfest01 = ln / time * reg
Else
Stehfest01 = 1
End If
End Function