Laplace 변환

라플라스 변환은 시간 평면 (t 평면)에서 신호를 s 평면으로 가져 오는 것으로,
미분 방정식을 대수적으로 풀기 위해 매우 효과적인 방법입니다.
예를 들면, 신호 처리에 있어서는 시간 변화하는 주기 신호를 라플라스 변환해, 다양한 조작(미적분 등 다수)을 행한 후에 역 라플라스 변환하여 시간축으로 되돌립니다.

심피

라플라스 변환을 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파의 라플라스 변환은 정의대로 계산해 이렇게 되고, 기억하고 있는 사람도 있다고 생각합니다.

삼각파

실제로 하고 싶은 것의 1보전에, 단순한 삼각파를 취급합니다.

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. 단위 신호를 생성(이 경우에는 삼각파 하나)
2. 라플라스 변환
3. 식 (1)을 곱한다
4. 역 라플러스 변환 수행

이 단계는 매우 쉬웠을 것입니다.
실제로 가보겠습니다.

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

자, 여기에서는 4개의 함수를 위에서부터 순서대로 만들었습니다.
첫 번째는 단위 신호의 라플라스 변환. (요소로 반환되므로 이후 Fs [0],
두 번째는 다음 식 (1)을 곱한 것입니다.
3번째는, 2번째를 적분하려고 한 것(s평면상에서는, 적분을 1/s를 곱하는 것으로 취급할 수 있다)
네 번째는 첫 번째를 적분하려고 시도한 것입니다.

중요한 것은, 특히 2번째, 3번째입니다.

그런데, 역변환해 보면,,,

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)"부분이 결과 출력을 방해하고 있습니다.
어떻게든 해결하고 싶습니다만, 경험이 있어 이해가 되는 분이 계시면 가르쳐 주었으면 합니다.