こんばんは、皆様、三頌亭です。数値的逆ラプラス変換のvbaによるユーザー定義関数のプログラムです。これは以前掲載しました有理関数の作例ですが、複素数計算の関数も併せて貼り付けるだけの形にしました。
Public Type Complex
x As Double
y As Double
End Type
Function ToComplex(x As Double, y As Double) As Complex
ToComplex.x = x
ToComplex.y = y
End Function
Function Cadd(z1 As Complex, z2 As Complex) As Complex
Cadd.x = z1.x + z2.x
Cadd.y = z1.y + z2.y
End Function
Function Csub(z1 As Complex, z2 As Complex) As Complex
Csub.x = z1.x - z2.x
Csub.y = z1.y - z2.y
End Function
Function Cmul(z1 As Complex, z2 As Complex) As Complex
Cmul.x = z1.x * z2.x - z1.y * z2.y
Cmul.y = z1.y * z2.x + z1.x * z2.y
End Function
Function Cdiv(z1 As Complex, z2 As Complex) As Complex
Dim c As Double
c = z2.x * z2.x + z2.y * z2.y
Cdiv.x = (z1.x * z2.x + z1.y * z2.y) / c
Cdiv.y = (z1.y * z2.x - z1.x * z2.y) / c
End Function
Function Conj(z As Complex) As Complex
Conj.x = z.x
comj.y = -z.y
End Function
Function Rez(z As Complex) As Double
Rez = z.x
End Function
Function Imz(z As Complex) As Double
Imz = z.y
End Function
Function Cabs(z As Complex) As Double
Cabs = Sqr(z.x * z.x + z.y * z.y)
End Function
Function atan2(x, y) As Double
Dim pi As Double
pi = 3.14159265358979
If x = 0# And y = 0# Then
atan2 = 0#
Exit Function
End If
If x = 0# Then
If y > 0 Then
atan2 = pi
Exit Function
Else
atan2 = -pi
Exit Function
End If
End If
If x > 0 Then
atan2 = Atn(y / x)
Else
If y > 0 Then
atan2 = pi + Atn(y / x)
Else
atan2 = -pi + Atn(y / x)
End If
End If
End Function
Function Cexp(z As Complex) As Complex
Cexp.x = Exp(z.x) * Cos(z.y)
Cexp.y = Exp(z.x) * Sin(z.y)
End Function
Function Clog(z As Complex) As Complex
Dim R As Double
R = z.x * z.x + z.y * z.y
If R <> 0 Then
Clog.x = 0.5 * Log(R)
Clog.y = atan2(z.x, z.y)
Else
Clog.x = 0#
Clog.y = 0#
End If
End Function
Function Cpow(z1 As Complex, z2 As Complex) As Complex
Cpow = Cexp(Cmul(Clog(z1), z2))
End Function
Function ccosh(z As Complex) As Complex
ccosh.x = (Exp(z.x) * Cos(z.y) + Exp(-z.x) * Cos(-z.y)) / 2
ccosh.y = (Exp(z.x) * Sin(z.y) + Exp(-z.x) * Sin(-z.y)) / 2
End Function
Function csinh(z As Complex) As Complex
csinh.x = (Exp(z.x) * Cos(z.y) - Exp(-z.x) * Cos(-z.y)) / 2
csinh.y = (Exp(z.x) * Sin(z.y) - Exp(-z.x) * Sin(-z.y)) / 2
End Function
Rem 数値的逆ラプラス変換filt07(t) Hosonoのアルゴリズム--------------
Rem atstp：分母 btstp：分子---ラプラス関数値の算出に必要なパラメータをセル参照
Rem 作例 F(s) = ht(8次多項式)/gt(8次多項式) : 有理関数
Public Function filt07(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(20, 0) As Double
Dim bred(20, 0) As Double
Dim pi As Double
Dim a01(20) As Complex
Dim b01(20) As Complex
Dim ht(200) As Complex
Dim gt(200) 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(500) As Double
Dim cp(500) 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 gt01:分母 ht01:分子
gt(7) = ss(n)
gt(6) = Cmul(gt(7), ss(n))
gt(5) = Cmul(gt(6), ss(n))
gt(4) = Cmul(gt(5), ss(n))
gt(3) = Cmul(gt(4), ss(n))
gt(2) = Cmul(gt(3), ss(n))
gt(1) = Cmul(gt(2), ss(n))
gt(0) = Cmul(gt(1), ss(n))
gt01 = a01(8)
gt01 = Cadd(gt01, Cmul(gt(7), a01(7)))
gt01 = Cadd(gt01, Cmul(gt(6), a01(6)))
gt01 = Cadd(gt01, Cmul(gt(5), a01(5)))
gt01 = Cadd(gt01, Cmul(gt(4), a01(4)))
gt01 = Cadd(gt01, Cmul(gt(3), a01(3)))
gt01 = Cadd(gt01, Cmul(gt(2), a01(2)))
gt01 = Cadd(gt01, Cmul(gt(1), a01(1)))
gt01 = Cadd(gt01, Cmul(gt(0), a01(0)))
ht(7) = ss(n)
ht(6) = Cmul(ht(7), ss(n))
ht(5) = Cmul(ht(6), ss(n))
ht(4) = Cmul(ht(5), ss(n))
ht(3) = Cmul(ht(4), ss(n))
ht(2) = Cmul(ht(3), ss(n))
ht(1) = Cmul(ht(2), ss(n))
ht(0) = Cmul(ht(1), ss(n))
ht01 = b01(8)
ht01 = Cadd(ht01, Cmul(ht(7), b01(7)))
ht01 = Cadd(ht01, Cmul(ht(6), b01(6)))
ht01 = Cadd(ht01, Cmul(ht(5), b01(5)))
ht01 = Cadd(ht01, Cmul(ht(4), b01(4)))
ht01 = Cadd(ht01, Cmul(ht(3), b01(3)))
ht01 = Cadd(ht01, Cmul(ht(2), b01(2)))
ht01 = Cadd(ht01, Cmul(ht(1), b01(1)))
ht01 = Cadd(ht01, Cmul(ht(0), b01(0)))
x1L(n) = Cdiv(ht01, gt01)
Rem ラプラス変換関数 Cadd Csub Cmul Cdiv ToComplex
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 変換値出力
filt07 = resigma * Exp(aa) / time
Else
Rem t=0の関数値：初期値の出力
filt07 = 0
End If
End Function
Rem 追加 平方根
Function csqr(z As Complex) As Complex
Dim xx As Double
Dim yy As Double
xx = z.x
yy = z.y
R = Sqr((xx) ^ 2 + (yy) ^ 2)
csqr.x = Sqr(xx + R) / Sqr(2)
csqr.y = Sqr(R - xx) / Sqr(2)
If z.y < 0 Then csqr.y = -csqr.y
End Function
