こんばんは、皆様、三頌亭です。今日は数値的な逆ラプラス変換をエクセルに関数として組み込むプログラムを紹介いたします。もともとこのプログラムは実験データを拡散方程式の解に当てはめるために作成したものです。もうだいぶ前に作成したものですが、最近になっても取っつきやすいエクセル上で使用できるものがないようなので公開いたします。とはいうもののそんな大げさな、かっこいいものではありません(^^;)。しかし実戦使用でだいぶ使ってきましたので、愛着のあるプログラムですw。
ところで複素数計算のコードは下記のhpよりお借りして、よく使う双曲線関数と平方根の関数を追加して作成してあります。
Fallen Physicist, Rising Engineer
https://sci.tea-nifty.com/blog/2008/10/excelvba-3b05.html?optimized=0
まずこのパッケージをエクセルから開発>visual basic>挿入>標準モジュールで貼り付けてください。そののちに逆変換のfilt06を追加で貼り付けてxlmsで保存して出来上がりです。
filt06(time As Double, atstp As Variant, btstp As Variant, nts As Integer) となっております。関数計算に必要なパラメータ・セットはatstpとbtstpでセル範囲を指定します。作例では有理関数の分母分子のsの8次式までの係数です。Ntsはオイラーの収束加速をもちない項数を入力します。20-40ぐらいでしょうか?。オイラーの収束加速は10項で固定してありますが、下記の係数のプログラムを組み込むことで大きくすることは可能です。
Public function eulercoeff(p as double, eff as double) as double
dim ap(100) as double
dim cp(100) as double
cp(0)=1
ap(p+1)=0
for q=1 to p
cp(q)=cp(q-1)*(p+1-q)/q
next q
for q=p to 0 step -1
ap(q)=ap(q+1)+cp(q)
next q
eulercoeff = ap(eff)
end function
数値的逆ラプラス変換のアルゴリズムはたくさんのものがあります。Stehfestのものをはじめとしていくつかのものを作成してみましたが、細野のものが最もロバストだったので現在もこれを使用しています。もちろん、時間と空間のメッシュを作成して実関数のみを使用する近似解法(ルンゲクッタ、クランク-ニコルソン)もあるのですが実験データの当てはめ(非線形最小二乗法)との相性があまりよくないのでこの方法に落ち着きました。ご参考までに追記しておきます。
参考文献
BASICによる高速ラプラス変換(細野敏夫:共立出版:1984)