こんばんは、皆様、三頌亭です。前回はあまり得意でない矩形波を逆変換した例をお見せいたしましたので、今回はベッセル関数を例にとって細野法で数値的逆ラプラス変換行ってエクセルの関数値と比較してみましょう。ところで最近ではフリーの数式処理システム（maximaなど）もありますので、それらで逆ラプラス変換をやってみて時間次元の解析解が簡単にに求まるならば数値計算にはそれを使えばよろしいです。逆変換が簡単な初等関数の組み合わせで表現できない場合、近似値を求める手段として意味があると三頌亭は考えます。
さて、ベッセル関数計算用のコードは以下の通りです。
Rem 数値的逆ラプラス変換filt22(t) Hosonoのアルゴリズム
Rem atstp：分母 btstp：分子---ラプラス関数値の算出に必要なパラメータをセル参照
Rem 作例 F(s) = (sqrt(s^2+1)-s)^n/sqrt(s^2+1) 第一種ベッセル関数Jn(t) :n=a01(0)
Public Function filt22(time As Double, atstp As Variant, btstp As Variant, nts As Integer, p As Double) As Double
Dim ss(10000) As Complex
Dim x1L(10000) As Complex
Dim res(10000) As Double
Dim aa As Double
Dim n As Double
Dim cnt As Double
Rem パラメータ関連・複素数変数
Dim resigma As Double
Dim ared(100, 1) As Double
Dim bred(100, 1) As Double
Dim pi As Double
Dim a01(20) As Complex
Dim b01(20) As Complex
Dim ht(20) As Complex
Dim gt(20) 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
Dim ap(1000) As Double
Dim cp(1000) As Double
Rem 各種定数
pi = 3.14159265358979
Rem オイラー変換係数
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
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 = 10
For n = 1 To nts + p
ss(n) = ToComplex(aa / time, (n - 0.5) / time * pi)
Rem ラプラス変換関数記述 Cadd Csub Cmul Cdiv ToComplex s = ss(n):
gt01 = csqr(Cadd(Cmul(ss(n), ss(n)), ToComplex(1, 0)))
ht01 = Cpow(Csub(gt01, ss(n)), a01(0))
x1L(n) = Cdiv(ht01, gt01)
Rem ラプラス変換関数記述 Cadd Csub Cmul Cdiv ToComplex s = ss(n)
res(n) = Imz(x1L(n)) * (-1) ^ n
Next
Rem nts項分加算
resigma = 0
For n = 1 To nts
resigma = resigma + res(n)
Next
Rem オイラー変換-p項分追加
For i = 1 To p
resigma = resigma + res(nts + i) * ap(i) / 2 ^ (p)
Next
Rem 変換値出力
filt22 = resigma * Exp(aa) / time
Else
Rem t=0の関数値：初期値の出力
filt22 = 0
End If
End Function
Rem https://kms130.hatenablog.com/entry/2020/12/03/232117
エクセルの内蔵の関数BESSELJでの計算と比較して、オイラー変換10項、未変換30項で相対誤差は10＾(-7)から10^(-8)程度です。０次から8次ぐらいまでやってみましたが非常に性能のいい近似式です。ｔ-＞∞できれいに収束するような関数に対しては素晴らしく近似精度がよろしいです。因みに↑のコードのaa は 8にしてあってだいたい有効桁7,8桁ということでほぼ注文通り収まってます（グラフでは重なってわかりません）。ラプラス変換した関数の記述部分は以下の3行のみです。
Rem F(s) = (sqrt(s^2+1)-s)^n/sqrt(s^2+1)： 第一種ベッセル関数Jn(t)
gt01 = csqr(Cadd(Cmul(ss(n), ss(n)), ToComplex(1, 0)))
ht01 = Cpow(Csub(gt01, ss(n)), a01(0))
x1L(n) = Cdiv(ht01, gt01)
