[PYTHON] La conversion de Sympy Laplace n'est pas pratique

Conversion de Laplace

La conversion de Laplace consiste à amener le signal dans le plan temporel (plan t) vers le plan s. C'est un moyen très efficace pour résoudre algébriquement les équations différentielles. Par exemple, dans le traitement du signal, un signal périodique à changement de temps est converti par Laplace, diverses opérations (telles que l'intégration fine, etc.) sont effectuées, puis une conversion de Laplace inverse est effectuée pour revenir à l'axe des temps.

Sympy Je voulais faire une conversion Laplace avec Python, et j'ai trouvé que cela pouvait être fait avec Sympy, alors je l'ai essayé, mais je ne pouvais pas l'utiliser à mort, donc je l'enregistrerai ici.

Connaissance de base de la conversion Laplace

Définition

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

Exemple simple

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

C'est un bon résultat. c'est correct. La transformation de Laplace de l'onde cos est calculée comme défini, et je pense que certaines personnes s'en souviennent.

Onde triangulaire

Une étape avant ce que vous voulez réellement faire, traitez une simple onde triangulaire.

TriangularPulse.py

T = 0.2 #période
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" est une fonction du côté lourd, qui est comme une fonction d'étape. Le signal ft a été tracé sur l'axe des temps et une onde triangulaire a été générée avec succès. La bonne chose à propos de la conversion de Laplace est que lorsque vous voulez la rendre continue périodiquement, sur le plan s

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

C'est un bon endroit pour postuler. En d'autres termes

1. 1. Générer un signal unitaire (dans ce cas, une onde triangulaire)
2. Conversion de Laplace
3. 3. Multiplier l'équation (1)
4. Effectuer une conversion de Laplace inverse

Cette étape aurait dû être très simple. J'irai vraiment.

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

Maintenant, ici, j'ai créé quatre fonctions dans l'ordre du haut. La première est la conversion de Laplace du signal unitaire. (Puisqu'il est retourné comme un élément, Fs [0], c'est-à-dire que la formule est extraite et utilisée à partir de maintenant) Le second est celui qui est ensuite multiplié par l'équation (1). Le troisième est une tentative d'intégration du second (sur le plan s, l'intégration peut être gérée en multipliant par 1 / s). Le quatrième est une tentative d'intégration du premier

Les deuxième et troisième sont particulièrement importants.

Au fait, quand j'essaye la conversion inverse, ...

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)

C'est devenu comme ça. Une fois tracés dans l'ordre, le premier et le quatrième ont été créés à ce stade.

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

est devenu. Le premier est, bien sûr, de revenir au signal d'origine. Le second semble être bon car le signal qui l'a intégré est sorti.

Partie d'erreur

Le deuxième signal est lorsque la conversion inverse est effectuée.

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)

Et la partie de "InverseLaplaceTransform (1 / (s ** 2 --s ** 2 * exp (-0.2 * s)), s, t, _None)" interfère avec la sortie du résultat. J'aimerais le résoudre d'une manière ou d'une autre, mais si quelqu'un a de l'expérience et comprend, faites-le moi savoir.