こんばんは、皆様、三頌亭です。細野法でオイラー変換の項数と未変換の項数を自由に設定できるようにしたVBAのコードです。これでおそらくラプラス変換した関数部分の記述を変えるだけでオッケーでしょうＷ。これらの数値的逆ラプラス変換法のアルゴリズムがあまり得意ではない方形波パルスの逆変換を試してみます。オイラー変換の項数と未変換の項数は20,20のものと80,400のもののグラフをお示しいたします。
Rem 数値的逆ラプラス変換filt10(t) Hosonoのアルゴリズム----------------
Rem atstp：分母 btstp：分子-ラプラス関数値の算出に必要なパラメータをセル参照
Rem 作例 F(s) = exp(-a*s)/s*(1-exp(-T*s)) 方形波パルス
Rem パラメータ入力：a01(0)=a,b01(0)=T
Public Function filt10(time As Double, atstp As Variant, btstp As Variant, nts As Integer, p As Double) As Double
Dim ss(5000) As Complex
Dim x1L(5000) As Complex
Dim res(5000) As Double
Dim aa As Double
Dim n As Integer
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(200) As Double
Dim cp(200) 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 = 8
For n = 1 To nts + p
ss(n) = ToComplex(aa / time, (n - 0.5) / time * pi)
Rem ラプラス変換関数記述 Cadd Csub Cmul Cdiv ToComplex F(s : ss(n))
ht01 = Cdiv(ToComplex(1,0),ss(n))
gt02 = Csub(ToComplex(0,0),a01(0))
gt02 = Cmul(ss(n),gt02)
gt01 = Csub(ToComplex(1,0),Cexp(gt02))
gt03 = Csub(ToComplex(0, 0), b01(0))
gt03 = Cmul(ss(n), gt03)
x1L(n) = Cmul(Cmul(ht01, gt01), Cexp(gt03))
Rem ラプラス変換関数記述 Cadd Csub Cmul Cdiv ToComplex F(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 変換値出力
filt10 = resigma * Exp(aa) / time
Else
Rem t=0の関数値：初期値の出力
filt10 = 0
End If
End Function
どうでしょう。まあ相変わらず無駄の多いみっともないコードで申し訳ないです（＾＾；）。ギブズ現象がはっきりわかります。時間を大きくしていくともっとはっきりしてきます。私の場合、拡散方程式でしたので急激な変動があまりありません。なので、あんまりわかりませんでした.。ルンゲクッタ法で近似解を求めても解はほとんど同じでした。因みによく論文なんかで使われているベッセル関数なんかは大体大丈夫です。ご参考までに・・・。
追記：時間大きくすると誤差が大きくなるという意味がわかりづらいと思いましたので、グラフの例を二つ追加いたしました。矩形波の連続パルスです。
作例 F(s) = (1/s)*(1/(1+exp(-T*s))) 方形波パルス（周期波）で関数部分のコードは
ht01 = Cdiv(ToComplex(1, 0), ss(n))
gt01 = Cadd(ToComplex(1, 0), Cexp(Cmul(ss(n), Csub(ToComplex(0, 0), a01(0)))))
gt01 = Cdiv(ToComplex(1, 0), gt01)
x1L(n) = Cmul(ht01, gt01)
です。
