[PYTHON] Sympy Laplace transform is not practical

Laplace transform

The Laplace transform brings the signal in the time plane (t plane) to the s plane. It is a very effective means for solving differential equations algebraically. For example, in signal processing, a time-changing periodic signal is Laplace-transformed, various operations (such as calculus, etc.) are performed, and then the inverse Laplace transform is performed to return to the time axis.

Sympy I wanted to do the Laplace transform with Python, and I found that it could be done with Sympy, so I tried it, but I couldn't use it to death, so I'll record it here.

Basic knowledge of Laplace transform

Definition

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

Simple example

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))

This is a good result. that's correct. The Laplace transform of the cos wave is calculated as defined, and I think some people remember it.

Triangle wave

One step before what you actually want to do, deal with a simple triangle wave.

TriangularPulse.py

T = 0.2 #period
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" is a heavyside function, which is like a step function. The signal ft was plotted on the time axis and a triangular wave was successfully generated. The good thing about the Laplace transform is that when you want to make it continuous periodically, on the s plane

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

It is a good place just to apply. In other words

1. 1. Generate a unit signal (in this case, one triangular wave)
2. Laplace transform
3. 3. Multiply equation (1)
4. Perform inverse Laplace transform

This step should have been very easy. I will actually go.

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

Now, here, I created four functions in order from the top. The first is the Laplace transform of the unit signal. (Since it is returned as an element, Fs [0], that is, the equation is extracted and used from now on) The second is the one that is then multiplied by equation (1). The third is an attempt to integrate the second (on the s plane, the integral can be handled by multiplying it by 1 / s). The fourth is an attempt to integrate the first

Of particular importance are the second and third.

By the way, when I try the reverse conversion, ...

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)

It has become like this. When plotted in order, the first and fourth were made at this stage.

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")

have become. The first is, of course, to return to the original signal. The second one seems to be good because the signal that integrates this came out.

Error part

The second signal is when the inverse conversion is performed.

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)

And the part of "InverseLaplaceTransform (1 / (s ** 2 --s ** 2 * exp (-0.2 * s)), s, t, _None)" is interfering with the output of the result. I'd like to solve it somehow, but if anyone has any experience and understands, please let me know.