SympyのLaplace変換が実用的じゃない

14905 ワード
Signalprocessing Python sympy Python テキストリンク

Laplace変換

ラプラス変換は，時間平面（t平面）における信号をs平面に持ってくることで，
微分方程式を代数的に解くために非常に有効な手段です．
例えば，信号処理においては時間変化する周期信号をラプラス変換して，色々な操作（微積分など他多数）を行ったのちに逆ラプラス変換して時間軸に戻します．

Sympy

ラプラス変換をPythonで行おうと思い，Sympyでそれができることがわかったので行ってみましたがこれが死ぬほど使えませんでしたのでここに記録しておきます．

ラプラス変換の基礎知識

定義


\mathcal{L}[f(t)] = \int_0^{\infty}f(t)e^{-st}dt

簡単な例

test.py
import sympy as sp

s, t = sp.symbols("s, t")
w = sp.symbols("w", real=True)
sp.laplace_transform(sp.cos(w*t), t, s)

>> (s/(s**2 + w**2), 0, Eq(Abs(periodic_argument(polar_lift(w)**2, oo)), 0))

これは真っ当な結果で良いです．正しいです．cos波のラプラス変換は定義通り計算してこうなりますし，覚えている人もいると思います．

三角波

実際に行いたいものの１歩前に、単純な三角波を扱います。

TriangularPulse.py

T = 0.2 #周期
ft = 2/T*t*sp.Heaviside(t) + 2/T*(t-T)*sp.Heaviside(t-T) - 4/T*(t-T/2)*sp.Heaviside(t-T/2)

p = plot(ft, (t, -0.2, 1), xlabel="t", ylabel="y")

「Heaviside」はヘヴィサイド関数のことで、ステップ関数みたいなもんです.
信号ftを時間軸でプロットし、うまく三角波を生成しました。
ラプラス変換の良いところは、これを周期的に連続させたいときにs平面上で


\frac{1}{1-e^{-Ts}}　...(1)

をかけるだけで良いところです。つまり

１．単位信号を生成（今の場合では三角波１つ）
２．ラプラス変換する
３．式(1)をかける
４．逆ラプラス変換を行う

このステップがとても楽なはずでした。
実際に行ってみます。

LaplaceTransform.py

ft = 2/T*t*sp.Heaviside(t) + 2/T*(t-T)*sp.Heaviside(t-T) - 4/T*(t-T/2)*sp.Heaviside(t-T/2)

Fs = laplace_transform(ft, t, s)
Fs_period = (1-sp.exp(-T*s))**(-1)*Fs[0]
Fs_period_integral = Fs_period / s
Fs_integral = Fs[0] / s
print(Fs)
print(Fs[0])
print(Fs_period)
print(Fs_period_integral)
print(Fs_integral)



(10.0/s**2 + 10.0*exp(-0.2*s)/s**2 - 20.0*exp(-0.1*s)/s**2, 0, True)
10.0/s**2 + 10.0*exp(-0.2*s)/s**2 - 20.0*exp(-0.1*s)/s**2
(10.0/s**2 + 10.0*exp(-0.2*s)/s**2 - 20.0*exp(-0.1*s)/s**2)/(1 - exp(-0.2*s))
(10.0/s**2 + 10.0*exp(-0.2*s)/s**2 - 20.0*exp(-0.1*s)/s**2)/(s*(1 - exp(-0.2*s)))
(10.0/s**2 + 10.0*exp(-0.2*s)/s**2 - 20.0*exp(-0.1*s)/s**2)/s

さて、ここでは、４つの関数を上から順につくりました。
１つ目は、単位信号のラプラス変換. (要素として返されるので、以後ではFs[0]、つまり式を取り出して使用した)
２つ目は、その後に式（１）をかけたもの
３つ目は、２つ目を積分しようとしたもの（s平面上では、積分を1/sをかけることで扱える）
４つ目は、１つ目を積分しようとしたもの

重要なのは、とくに２つ目、３つ目です。

さて、逆変換してみると、、、

InverseLaplaceTransform.py

ft = sp.inverse_laplace_transform(Fs[0],s,t,noconds=True)
ft_period = sp.inverse_laplace_transform(Fs_period,s,t)
ft_period_integral = sp.inverse_laplace_transform(Fs_period_integral,s,t)
ft_integral = sp.inverse_laplace_transform(Fs_integral,s,t)
print(ft)
print("******************")
print(ft_period)
print("******************")
print(ft_period_integral)
print("******************")
print(ft_integral)


10.0*t*Heaviside(t) + (10.0*t - 2.0)*Heaviside(t - 0.2) - (20.0*t - 2.0)*Heaviside(t - 0.0999999999999999)
******************
10.0*InverseLaplaceTransform(1/(s**2 - s**2*exp(-0.2*s)), s, t, _None) - 20.0*InverseLaplaceTransform(1/(s**2*exp(0.1*s) - s**2*exp(-0.1*s)), s, t, _None) + 10.0*InverseLaplaceTransform(1/(s**2*exp(0.2*s) - s**2*exp(-5.55111512312578e-17*s)), s, t, _None)
******************
10.0*InverseLaplaceTransform(1/(s**3 - s**3*exp(-0.2*s)), s, t, _None) - 20.0*InverseLaplaceTransform(1/(s**3*exp(0.1*s) - s**3*exp(-0.1*s)), s, t, _None) + 10.0*InverseLaplaceTransform(1/(s**3*exp(0.2*s) - s**3*exp(-5.55111512312578e-17*s)), s, t, _None)
******************
5.0*t**2*Heaviside(t) + (5.0*t**2 - 2.0*t + 0.2)*Heaviside(t - 0.2) - (10.0*t**2 - 2.0*t + 0.0999999999999998)*Heaviside(t - 0.0999999999999999)

こうなってしまいました。順にプロットしたところ、この段階でできたのは１つ目と４つ目で

Plot1.py
p1 = plot(ft, (t, -0.2, 1), xlabel="t", ylabel="y")
p4 = plot(ft_integral, (t, -0.2, 1), xlabel="t", ylabel="y")


となりました。１つ目は、もとの信号に戻って当然です。２つ目はこれを積分した信号が出てきたので良さそうです。

エラー部分

２つ目の信号は、逆変換を行った際に

10.0*InverseLaplaceTransform(1/(s**2 - s**2*exp(-0.2*s)), s, t, _None) - 20.0*InverseLaplaceTransform(1/(s**2*exp(0.1*s) - s**2*exp(-0.1*s)), s, t, _None) + 10.0*InverseLaplaceTransform(1/(s**2*exp(0.2*s) - s**2*exp(-5.55111512312578e-17*s)), s, t, _None)

となってしまい、"InverseLaplaceTransform(1/(s*2 - s*2*exp(-0.2*s)), s, t, _None)"の部分が結果の出力を邪魔しています。
どうにか解決したいのですが、経験があっておわかりになる方がいらっしゃいましたら教えてほしいです。