数値的逆ラプラス変換(Gaver-Stehfest 法について:vba-code作例)

こんばんは、皆様、三頌亭です。数値的逆ラプラス変換の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

 

f:id:kms130:20210114212603p:plain

f(t)=exp(-alpha*t) : F(s)=1/(s+alpha)

 

f:id:kms130:20210114224738j:plain

Stehfest table01

f:id:kms130:20210114224828j:plain

Stehfest table02