数値的逆ラプラス変換（Gaver-Stehfest 法について）

こんばんは、皆様、三頌亭です。このまえから解説してます逆ラプラス変換方法ですが、ご参考までに細野法ではない、Gaver-Stehfest 法について簡単にご説明いたします。まあ大体、利用者の立場としては「能書きはいいから早くプログラムをよこせ！」というのがよくあるお話ですねｗｗ。

ところでGaver-Stehfest 法についてですが複素数を必要とせず、実数のみで計算が完了いたします。ラプラス変換済みの関数F(s)にそのまま実数s=Ln(2)*n/time(n=1～20程度まで)を入力してF(s)の関数値を計算いたします。それぞれの値に定まった係数をかけて加算していきます。総和にLn(2)/timeを掛けて、timeにおける時間次元の関数値f(time)とします。因みにかける係数は下記のとおりです。一番下の画像に一般式をお示します。

vv(1) = -5.51146384479718E-06

vv(2) = 0.152386463844797

vv(3) = -117.465476190476

vv(4) = 17342.4493386243

vv(5) = -922806.928902116

vv(6) = 23774087.7871031

vv(7) = -349421166.19537

vv(8) = 3241369852.23187

vv(9) = -20276948307.2378

vv(10) = 89464829823.7973

vv(11) = -287020921147.103

vv(12) = 682992010281.513

vv(13) = -1219082330054.38

vv(14) = 1637573800842.02

vv(15) = -1647177486836.12

vv(16) = 1221924554444.23

vv(17) = -648806558817.535

vv(18) = 233316653213.707

vv(19) = -50913800705.4676

vv(20) = 5091380070.54676

値をご覧になればお分かりになると思いますが、振動してどんどん大きくなりますｗ。なのでｎを大きくとると相対誤差が非常に悪くなってしまいます。一応、一般のｎについて係数を計算するvbaのコードを紹介いたします（金融工学の論文からのコピーです。貼り付ければそのままで動作します）。

https://www.researchgate.net/publication/228310935_Option_Pricing_with_Jumps

一番最後の部分がそうです。

Rem StehfestCoefficients returns coefficients in Stehfest inversion formula
Function StehfestCoefficients(Nterms As Integer) As Variant
Dim i As Integer, j As Integer, k As Integer, nh As Integer, sn As Integer, nn As Integer
nh = Nterms / 2 'Must be even
Nterms = 2 * nh
ReDim Coeffs(0 To Nterms) As Double, g(0 To Nterms) As Double, h(0 To nh) As Double
g(0) = 1
For i = 1 To Nterms
g(i) = g(i - 1) * i
Next i
h(1) = 2# / g(nh - 1)
For i = 2 To nh
h(i) = Exp(nh * Log(i)) * g(2 * i) / (g(nh - i) * g(i) * g(i - 1))
Next i
If nh Mod 2 = 0 Then sn = 1 Else sn = -1
For i = 1 To Nterms
Coeffs(i) = 0
If i < nh Then nn = i Else nn = nh
j = Int( (i + 1) / 2)
For k = j To nn
Coeffs(i) = Coeffs(i) + h(k) / (g(i - k) * g(2 * k - i))
Next k
sn = -sn
Coeffs(i) = Coeffs(i) * sn
Next i
StehfestCoefficients = Coeffs
End Function