4. Linearni časovno invariantni sistemi z eno prostostno stopnjo#

Sistema z eno prostostno stopnjo (ang. single degree of freedom system - SDOF, slika spodaj) lahko opišemo z nehomogeno navadno diferencialno enačbo drugega reda:

\[ m\,\ddot{x}(t)+c\,\dot{x}(t)+k\,x(t)=f(t), \]

structural_SDOF_model

kjer je \(f(t)\) vzbujevalna sila. Rešitev nehomogene enačbe je v obliki:

\[ x(t) = x_{h}(t) + x_{p}(t), \]

kjer je \(x_{h}\) homogena rešitev in \(x_{p}\) partikularna rešitev.

Prosti odziv (homogena rešitev)#

Če vzbujevalne sile ni, \(f(t)=0\), je homogena rešitev pridobljena iz:

\[ m\,\ddot{x}_{h}(t)+c\,\dot{x}_{h}(t)+k\,x_{h}(t)=0. \]

Uporabimo nastavek za rešitev \(x_{h}(t)\):

\[ x_{h}(t)=X_{h}\,\mathrm{e}^{\lambda \,t}, \]

kjer sta amplituda premika \(X_{h}\) in lastna vrednost \(\lambda\) neznani (in ju je treba določiti). Vstavimo nastavek v diferencialno enačbo:

\[ \lambda^2\,m\,X_{h}\,\mathrm{e}^{\lambda \,t}+\lambda\,c\,X_{h}\,\mathrm{e}^{\lambda \,t}+k\,X_{h}\,\mathrm{e}^{\lambda \,t}=0. \]

Izraz uredimo:

\[ \left(\lambda^2\,m+\lambda\,c+k\right)\,X_{h}\,\mathrm{e}^{\lambda \,t}=0. \]

Netrivialno (\(X_{h}\,\mathrm{e}^{\lambda \,t}\neq0\)) rešitev določimo na podlagi karakteristične enačbe:

\[ \lambda^2\,m+\lambda\,c+k=0. \]

Rešitev karakteristične enačbe je:

\[ \lambda_{1,2}=-\frac{c}{2\,m} \pm \frac{\sqrt{c^2-4\,k\,m}}{2\,m}, \]

kjer sta \(\lambda_{1,2}\) lastni vrednosti sistema (v splošnem kompleksno število).

Glede na razmerje med vztrajnostno, elastično in dušilno silo je rešitev:

  • nadkritično dušeno: če sile dušenja prevladajo nad vztrajnostnimi in prožnostnimi sile (\(c^2\geq 4\,k\,m\)),

  • kritično dušeno: če so sile dušenja enake vztrajnostnim in prožnostnim silam (\(c^2=4\,k\,m\)),

  • podkritično dušeno: če vztrajnostne in prožnostne sile prevladajo nad dušilnimi (\(c^2\leq 4\,k\,m\)).

V primeru nadkritično in kritično dušenega sistema je rešitev v realni obliki in gibanje ni oscilatorno.

V okviru te knjige se osredotočamo na podkritično dušen primer, ko so kompleksne lastne vrednosti konjugirani pari. Masa \(m\), dušenje \(c\) in togost \(k\) pa definirajo nedušeno lastno krožno frekvenco (tudi resonančna frekvenca, lastna frekvenca):

\[ \omega_0 = \sqrt{\frac{k}{m}}. \]

Tukaj je pomembno spomniti (glejte osnovne predmete na temo dinamike in mehanskih nihanj), da pri linearnih sistemih prednapetje ne spremeni lastne frekvence.

Dušenje je mogoče definirati s pomočjo razmernika dušenja:

\[ \delta = \frac{c}{c_{c}}, \]

kjer je:

\[ c_{c} = 2\,\sqrt{k\,m} \]

kritično dušenje.

S pomočjo nedušene lastne frekvence \(\omega_0\) in razmernika dušenja \(\delta\) lahko homogeno diferencialno enačbo preoblikujemo v standardno obliko:

\[ \ddot{x}_{h}(t)+2\,\delta\,\omega_0\,\dot{x}_{h}(t)+\omega_0^2\,x_{h}(t) =0. \]

Nadalje \(\omega_0\) in \(\delta\) vstavimo v izraz za lastne vrednosti:

\[ \lambda_{1,2}=-\delta\, \omega_0 \pm \mathrm{i}\,\omega_{d}, \]

kjer se:

\[ \omega_{d} = \omega_0\,\sqrt{1-\delta^2} \]

imenuje dušena lastna krožna frekvenca (pravilno jo označimo z \(\omega_{0,d}\), vendar tukaj zaradi zapisa v kodi, spodaj, uporabimo kratko obliko).

Izpeljane lastne vrednosti vstavimo v nastavek za homogeno rešitev \(x_{h}(t)\). Ker imamo dve lastni vrednosti, je splošna rešitev linearna kombinacija obeh:

\[\begin{split} \begin{split} x_{h}(t)&=X_{h,1}\,\mathrm{e}^{\lambda_1 \,t}+X_{h,2}\,\mathrm{e}^{\lambda_2 \,t}\\ &=\mathrm{e}^{-\delta\, \omega_0\,t}\, \big( X_{h,1}\,\mathrm{e}^{+ \mathrm{i}\, \omega_{d}\,t}+X_{h,2}\,\mathrm{e}^{- \mathrm{i}\, \omega_ {d} \,t} \big), \end{split} \end{split}\]

kjer sta konstanti \(X_{h,1}\) in \(X_{h,2}\) odvisni od začetnih pogojev.

Pri impulznem odzivu bomo gradili na začetnih pogojih, ko sistem izpustimo iz ravnotežja s hitrostjo različno od nič:

\[\begin{split} \begin{aligned} x(0)&=0\to& X_{h,1}+X_{h,2}&=0\\ \dot x(0)&=\dot x_0\to& X_{h,1}\,\lambda_1+X_{h,2}\,\lambda_2&=\dot x_0. \end{aligned} \end{split}\]

Rešitev začetnih pogojev predstavlja t.i. impulzno odzivno funkcijo.

Impulzna prenosna funkcija#

Opomba

Impulzna prenosna funkcija:

\[ x_{h}(t)=\frac{\dot x_0}{\omega_{d}}\,\mathrm{e}^{-\delta\,\omega_0\,t}\,\sin(\omega_{d}\,t). \]

Odziv je sestavljen iz ovojnice \({\dot x_0}/{\omega_{d}}\,\exp(-\delta\, \omega_0\,t)\) in nihajočega dela \(\sin(\omega_{d}\,t)\), glejte sliko spodaj.

Hide code cell source
import sympy as sym
t, ω_0, ω_d, d_x_0, δ= sym.symbols('t, \\omega_0, \\omega_d, d_x_0, \\delta', real=True)
podatki = {ω_0: 1, d_x_0: 1, δ: 0.1}

ω_d = ω_0 * sym.sqrt(1-δ**2)

ovojnica = d_x_0/ω_d * sym.exp(-δ*ω_0*t)

x = ovojnica * sym.sin(ω_d*t)

p1 = sym.plot(ovojnica.subs(podatki), (t,0,50), line_color='C0', xlabel='$t$ [s]', ylabel='$x_h(t)$', title='Časovna domena',
             label='Ovojnica', show=False)
p2 = sym.plot(x.subs(podatki), (t,0,50), line_color='C1', 
             label='$x_h(t)$', show=False)
p1.extend(p2)
p1.legend = True
p1.show()
../_images/994928dbb4f730971f057b6d36e4be7dc2cb90fafb5044adffd5f9bedaacdd46.png

Harmonsko vzbujanje#

Nadaljujemo z rešitvijo gibalne enačbe za primer harmonskega vzbujanja:

\[ f(t)=F\,\mathrm{e}^{ \mathrm{i}\,\omega\,t}=F\,\underbrace{\big(\cos(\omega\,t)+ \mathrm{i}\,\sin(\omega\,t)\big)}_{\mathrm{e}^{ \mathrm{i}\,\omega\,t}}, \]

kjer je \(F\) amplituda sile (v splošnem kompleksna vrednost) in \(\omega\) krožna frekvenca vzbujanja.

Partikulano rešitev \(x_p(t)\) določimo iz enačbe:

\[ m\,\ddot{x}_{p}(t)+c\,\dot{x}_{p}(t)+k\,x_{p}(t)=f(t), \]

tako, da predpostavimo harmonski odziv:

\[ x_{p}(t)=X_{p}\,\,\mathrm{e}^{ \mathrm{i}\,\omega\,t}, \]

kjer je \(X_{p}\) neznana kompleksna amplituda partikularne komponente (včasih imenovana tudi fazor). Nadaljujemo z izpeljavo:

\[ -\omega^2\,m\,X_{p}\,\mathrm{e}^{ \mathrm{i}\,\omega\,t}+ \mathrm{i}\,\omega\,c\,X_{p}\,\mathrm{e}^{ \mathrm{i}\,\omega\,t}+k\,X_{p}\,\,\mathrm{e}^{ \mathrm{i}\,\omega\,t}=F\,\mathrm{e}^{ \mathrm{i}\,\omega\,t}. \]

Netrivialna rešitev zahteva \(\mathrm{e}^{ \mathrm{i}\,\omega\,t}\neq 0\) in zato rešitev iščemo na podlagi enačbe:

\[ -\omega^2\,m\,X_{p}+ \mathrm{i}\,\omega\,c\,X_{p}+k\,X_{p}=F, \]

kar vodi v določitev razmerja med odzivom \(X_p\) in vzbujanjem \(F\), kar imenujemo tudi frekvenčna prenosna funkcija \(\alpha(\omega)\):

\[ \alpha(\omega)=\frac{X_{p}}{F}=\frac{1}{-\omega^2\,m+ \mathrm{i}\,\omega\,c+k}. \]

Če izpostavimo \(1/k\), nato uporabimo \(\omega_0^2=k/m\), \(c=\delta\,c_c\), \(c_c=2\sqrt{k\,m}\) ter dodatno vpeljemo razmernik frekvenc \(r=\omega/\omega_0\), izpeljemo:

\[ \alpha(r)=\frac{1}{k}=\frac{1}{1-r^2+2\,\mathrm{i}\,\delta\,r}. \]

V nadaljevanju izpeljani izraz prikažemo grafično.

Hide code cell source
import sympy as sym
r, δ, k= sym.symbols('r, \\delta, k', real=True)
i = sym.I

α = 1/k/(1-r**2 + i*2*δ*r)
podatki = {k: 1, δ:0.1}
podatki2 = {k: 1, δ:0.5}

p1 = sym.plot(sym.Abs(α.subs(podatki)), (r,0,3), line_color='C0', xlabel='$r$ [\\cdot]', ylabel='$|\\alpha(r)|$', 
             title=f'Amplituda frekvenčne prenosne funkcije',
             label=f'$\\delta=${δ.subs(podatki):2.1f}', show=False)
p2 = sym.plot(sym.Abs(α.subs(podatki2)), (r,0,3), line_color='C1', 
             label=f'$\\delta=${δ.subs(podatki2):2.1f}', show=False)
p1.extend(p2)
p1.legend = True
p1.show()

p1 = sym.plot(sym.arg(α.subs(podatki)), (r,0,3), line_color='C0', xlabel='$r$ [\\cdot]', ylabel='$\\angle\\alpha(\\omega)$', 
             title=f'Faza frekvenčne prenosne funkcije',
             label=f'$\\delta=${δ.subs(podatki):2.1f}', show=False)
p2 = sym.plot(sym.arg(α.subs(podatki2)), (r,0,3), line_color='C1', 
             label=f'$\\delta=${δ.subs(podatki2):2.1f}', show=False)
p1.extend(p2)
p1.legend = True
p1.show()
../_images/eb906278ddb7ab2ef839d0758caa3d383e8063510e33fddcd2bcf9c1d8fb0286.png ../_images/c80f982f108749a3779fc2b2b5cade0d74f71541da306af4e8375effaffa18e1.png

Partikularna rešitev je torej:

\[ x_{p}(t)=\frac{1}{-\omega^2\,m+ \mathrm{i}\,\omega\,c+k}\,F% \,\mathrm{e}^{ \mathrm{i}\,\omega\,t}. \]

Prehodni odziv in odziv v ustaljenem stanju.#

Najprej smo določili homogeno rešitev \(x_{h}(t)\), katera predstavlja prehodni odziv na neko začetno motnjo. Partikularna rešitev \(x_{p}(t)\) predstavlja rešitev na stacionarno harmonsko motnjo. Splošna rešitev je določena z vsoto obeh (slika spodaj).

sdof_homo_part_trans_steady

Frekvečna prenosna funkcija#

Najprej se bomo osredotočili zgolj na partikularno rešitev (vendar ne bomo več uporabljali podpisa \(p\)):

\[ \frac{X}{F}=\frac{1}{-\omega^2\,m+ \mathrm{i}\,\omega\,c+k}. \]

Ker zgornji izraz predstavlja frekvenčni odziv, se imenuje frekvenčna prenosna funkcija (ang. Frequency Response Function, tudi FRF, glejte Povezava med frekvenčno in impulzno prenosno funkcijo za opis načina, kako se FRF dobi iz homogene rešitve:

Opomba

Frekvenčna prenosna funkcija:

\[ \alpha(\omega)=\frac{1}{-\omega^2\,m+ \mathrm{i}\,\omega\,c+k}. \]

Z uporabo nedušene laste krožne frekvence \(\omega_0^2=k/m\), razmernika dušenja \(\delta=c/(2\,\sqrt{k\,m})\) in razmernika frekvenc ali tudi relativne frekvence \(r=\omega/\omega_0\), dobimo frekvenčni odziv v standardni obliki:

\[ \alpha(\omega)=\frac{1}{k}\cdot\frac{1}{1-r^2+2\, \mathrm{i}\,\delta\,r}. \]

\(\alpha(\omega)\) se imenuje frekvenčna prenosna funkcija ali bolj natančno dimenzijska frekvenčna prenosna funkcija (ang. dimensional frequency response function). V nasprotju z \(\alpha(\omega)\) je \(H_0(\omega)\) opredeljena kot:

\[ \begin{split} H_0(\omega)&=\frac{X(\omega)}{X_0} \end{split} \]

in je v nekaterih knjigah (npr.: Géradin and Rixen [2015]) imenovana brezdimenzijska frekvenčna prenosna funkcija (ang. nondimensional frequency response function). V izrazu zgoraj je statični poves \(X_0\) zaradi statične sile \(F\) definiran kot:

\[ X_0=\frac{F}{k}. \]

Opomba

Harmonski odziv \(X(\omega)\) lahko z brezdimenzionalno frekvenčno prenosno funkcijo definiramo kot:

\[\begin{split} \begin{split} H_0(\omega)&=\frac{X(\omega)}{X_0}\\ &=\frac{k}{-\omega^2\,m+ \mathrm{i}\,\omega\,c+k}\\ &=\frac{1}{1-r^2+ \mathrm{i}\,2\,\delta\,r}, \end{split} \end{split}\]

Oziroma:

\[ \alpha(\omega)=\frac{H_0(\omega)}{k}. \]

S pomočjo frekvenčne prenosne funkcije je mogoče odziv stacionarnega stanja pri harmonskem vzbujanju s \(\omega_f\) opisati kot:

\[ x(t)=\alpha(\omega=\omega_f)\,F\,\mathrm{e}^{ \mathrm{i}\,\omega_f\,t}. \]

Različne oblike FRF#

Dimenzijski FRF je smiseln, saj je statični odklon \(X_0\) težko izmeriti in se pogosto merita vzbujevalna sila in odzivni pospešek. Namesto pospeška pogosto merimo še hitrost in pomik. Posledično obstajajo različne oblike frekvenčne prenosne funkcije:

Opomba

\[\begin{split} \begin{aligned} \alpha(\omega)&=\frac{X(\omega)}{F(\omega)},&&\textrm{podajnost ali receptanca (ang. receptance),}\\ Y(\omega)&=\frac{\dot X(\omega)}{F(\omega)}= \mathrm{i}\,\omega\,\alpha(\omega),&&\textrm{pomičnost (ang. mobility),}\\ A(\omega)&=\frac{\ddot X(\omega)}{F(\omega)}=-\omega^2\,\alpha(\omega),&&\textrm{pospešenost (ang. accelerance),}\\ \end{aligned} \end{split}\]

kjer \(\dot X(\omega)\) in \(\ddot X(\omega)\) označujejo kompleksni amplitudi hitrosti ali pospeška v frekvenčni domeni. Odziv različih oblik FRF je prikazan na sliki spodaj:

  • podajnost: podresonančno konvergira k statičnemu inverzu togosti \(1/k\), nadresonančno pa kvadratično konvergira proti nič,

  • pomičnost: podresonančno linearno konvergira k nič, nadresonančno linearno konvergira proti nič,

  • pospešenost: podresonančno kvadratično konvergira proti nič, nadresonančno konvergira k inverzu mase \(1/m\).

Hide code cell source
import sympy as sym
ω, m, c, k= sym.symbols('\\omega, m, c, k', real=True)
i = sym.I

α = 1/(-ω**2 * m + i*ω*c + k)
podatki = {m: 1, c: 0.1, k: 1}

p1 = sym.plot(sym.Abs(α.subs(podatki)), (ω,0.01,100), line_color='C0', xlabel='$\\omega$ [rad/s]', ylabel='',  
             title=f'Amplituda frekvenčne prenosne funkcije ($\\omega_0=${sym.sqrt(k/m).evalf(subs=podatki):3.2f})',
             label='$|\\alpha(\\omega)|$ - podajnost', xscale='log', yscale='log', show=False)
p2 = sym.plot(sym.Abs(i*ω*α.subs(podatki)), (ω,0.01,100), line_color='C1', 
             label='$|Y(\\omega)|$ - pomičnost', xscale='log', yscale='log', show=False)
p3 = sym.plot(sym.Abs(-ω**2*α.subs(podatki)), (ω,0.01,100), line_color='C2', 
             label='$|A(\\omega)|$ - pospešenost', xscale='log', yscale='log', show=False)
p1.extend(p2)
p1.extend(p3)
p1.legend = True
p1.show()
../_images/77ba461cd9da03bfb526f6353b95c2039ce4d6eb6172711c0142ec41b69c9431.png

Opomba

Inverzna vrednost pospešenosti je dinamična masa (\(F/\ddot{X}\)), inverzna vrednost pomičnosti je mehanska impedanca (\( F/\dot{X}\)) in inverzna vrednost podajnosti je dinamična togost (\(F/X\)).

Lastnosti linearnega časovno invariantnega sistema#

Diferencialna enačba za sistem z eno prostostno stopnjo poveže vzbujanje (input) z odzivom (output):

\[ \underbrace{f(t)}_{\textrm{Vzbujanje}}\to {\fbox{Sistem}}\to \underbrace{x(t)}_{\textrm{Odziv}}. \]

Sistem je linearen, če ima lastnosti aditivnosti:

\[ f_1(t)+f_2(t) \to {\fbox{Linearen sistem}} \to x_1(t)+x_2(t) \]

in skalirnosti (tudi homogenost):

\[ a\,f(t) \to {\fbox{Linearen sistem}} \to a\,x(t) \]

Če je sistem linearen, lahko uporabimo načelo superpozicije:

\[ a_1\,f_1(t)+a_2\,f_2(t) \to {\fbox{Linearen sistem}} \to a_1\,x_1(t)+a_2\,x_2(t) \]

Še ena pomembna lastnost je časovna invariantnost:

\[ f(t-t_0) \to {\fbox{Časovno invarianten sistem}} \to x(t-t_0) \]

Odziv na periodično vzbujanje#

Zapis vzbujanja s Fourierovimi vrstami#

Glede Fourierove vrste v eksponentni (kompleksni) obliki lahko periodično vzbujanje \(f(t)\) zapišemo kot:

\[ f(t) = \sum_{n=-\infty}^{+\infty} c_n\,\mathrm{e}^{ \mathrm{i}\,n\,\omega_{p}\,t}, \]

kjer je \(\omega_{p}\) osnovna krožna frekvenca periodičnega vzbujanja:

\[ \omega_{p} = \frac{2\,\pi}{T}. \]

in so Fourierovi (kompleksni) koeficienti \(c_n\) za \(n\in \mathbb{Z}\) definirani kot:

\[ c_n = \frac{1}{T}\,\int_{\tau}^{\tau+T}f(t)\,\mathrm{e}^{- \mathrm{i}\,n\,\omega_{p}\,t}\,\textrm{d} t, \]

kjer je \(c_0\) povprečna vrednost \(f(t)\). Velja tudi \(c_{-n}=c_n^*\), kjer \({}^*\) označuje kompleksno konjugacijo.

Opomba

Ker velja princip superpozicije, lahko vzbujanje razstavimo na poljubno komponent; tukaj vzbujanje razstavimo na posamezne člene Fourierove vrste; Fourierov koeficient \(c_n\) tako definira eno harmonsko komponento.

Zapis odziva s pomočjo Fourierovih vrst#

Tukaj nadaljujemo izračunom rešitve za vsako harmonsko komponento vzbujanja posebej (glejte Frekvečna prenosna funkcija)

\[\begin{split} \begin{split} X_n&=X_{0,n}\,H_0(\omega_{p}\,n)\\ &=\frac{c_n}{k}\,H_0(\omega_{p}\,n)\\ &=c_n\,\alpha(\omega_{p}\,n), \end{split} \end{split}\]

kjer je \(c_n\) definiran zgoraj pri periodičnem vzbujanju s silo.

Ko imamo enkrat \(X_n\) določen, lahko izračunamo odziv te harmonske komponente v času:

\[ x_n(t)=X_n\,\,\mathrm{e}^{ \mathrm{i}\,n\,\omega_{p}\,t} \]

in se vrnemo k principu superpozicije ter odzive seštejemo:

\[ x(t) = \sum_{n=-\infty}^{+\infty} x_n(t). \]

Posplošitev periodičnega vzbujanja na poljubno vzbujanje#

Fourierova transformacija predstavlja posplošitev Fourierovih vrst tudi na neperiodične podatke. Izpeljave, ki smo jih s pomočjo Fourierovih vrst naredili za periodično vzbujanje, lahko tako posplošimo.

Opomba

Fourierova transformacija vzbujanja za prehod v frekvenčno domeno:

\[ F(\omega)=\mathcal{F}\{f(t)\} \]

vzbujanje razstavi na harmonske komponente \(F(\omega)\), ki jih nato rešimo po posameznih krožnih frekvencah \(\omega\):

\[ X(\omega)=\alpha(\omega)\,F(\omega) \]

in nato naredimo spet nazaj prehod v časovno domeno:

\[ x(t)=\mathcal{F}^{-1}\{X(\omega)\}. \]

Vzbujanje z enotskim impulzom#

Obravnavajmo sedaj vzbujanje enomasnega sistema z enotskim impulzom \(f_1(t)=\delta(t-\tau)\), kjer je \(\delta(t)\) Diracova delta funkcija (glejte: Diracova delta funkcija):

\[ \Delta I_1=\lim_{\Delta \tau\to 0}\int_{\tau}^{\tau+\Delta \tau}f_1(t)\,\textrm{d}t =1. \]

Enomasni sistem je ob času \(\tau^-\) v mirovanju, infinitezimalni čas pozneje \(\tau^+\) pa ima zaradi enotskega impulza \(\Delta I_1\) hitrost:

\[ \dot x_0 = \frac{\Delta I_1}{m}= \frac{1}{m}. \]

Po enotskem impulzu \(t\ge \tau\) je sila \(f_1(t)\) nič in sistem se odziva glede na homogeno rešitev pri ničelnem začetnem premiku \(x=0\) in \(\dot x=\dot x_0\) (glejte Impulzna prenosna funkcija):

\[ h(t-\tau)=x(t)=\frac{1}{m\,\omega_{d}}\,\mathrm{e}^{-\zeta\,\omega_0\,(t-\tau)}\,\sin(\omega_{d}\,(t-\tau)),\quad \textrm{za: } t>\tau, \]

kjer je \(h(t-\tau)\) impulzna prenosna funkcija.

Odziv na enotski impulz vzročnega sistema (tudi kavzalni sistem) je odvisen samo od preteklega (in trenutnega) vzbujanja, ne pa od prihodnjega. Vsi sistemi v tej knjigi veljajo za vzročne.

Posplošitev enotskega vzbujanja na poljubno vzbujanje#

Splošno vzbujevalno silo (slika spodaj) lahko vidimo kot superpozicijo enotskih impulzov, pomnoženih z določeno konstanto.

general_forcing

Vsak neskončno kratek impulz bi povzročil odziv glede na impulzno prenosno funkcijo \(h(t-\tau)\) pomnožen z amplitudo impulza \(f(\tau)\,\Delta\tau\). Za linearne sisteme take odzive po principu superpozicioniran seštejemo:

\[ x(t) = \sum_{\tau}f(\tau)\,h(t-\tau)\,\Delta \tau,\quad \textrm{za: } t>\tau, \]

zaradi infinitezimalno majhnih časovnih korakov \(\Delta \tau\to0\) vsota preide v integral:

\[ x(t) = \int_0^t\,f(\tau)\,h(t-\tau)\,\textrm{d} \tau,\quad \textrm{za: } t>\tau. \]

Zgornji izraz imenujemo konvolucijski integral, superpozicijski integral ali konvolucija (kdaj tudi Duhamelov integral).

Ker je impulzna prenosna funkcija \(h(t-\tau)=0\) za \(t<\tau\), se lahko spodnja meja integrala razširi v \(-\infty\):

\[ x(t) = \int_{-\infty}^t\,f(\tau)\,h(t-\tau)\,\textrm{d} \tau,\quad \textrm{za: } t>\tau. \]

Nadalje, če spremenimo spremenljivko integracije v \(\tau_1=t-\tau\) (posledično velja: \(\textrm{d} \tau_1=-\textrm{d} \tau\)), se meje integracije spremenijo v: \(-\infty\to+ \infty\) in \(t\to 0\):

\[ x(t) = -\int_{+\infty}^0\,f(t-\tau_1)\,h(\tau_1)\,\textrm{d} \tau_1 = \int_0^{+\infty}\,f(t-\tau_1)\,h(\tau_1)\,\textrm{d} \tau_1, \]

Sledi, da lahko povzamemo lastnost komutativnost konvolucije:

\[ x(t) = \int_{-\infty}^{t}\,f(\tau)\,h(t-\tau)\,\textrm{d} \tau = \int_{0}^{+\infty}\,h(\tau)\,f(t-\tau)\,\textrm{d} \tau, \]

ali v kratki obliki:

\[ x(t) = f(t) * h(t)= h(t)*f(t), \]

kjer \(*\) označuje konvolucijo. Primer (spodaj) prikazuje vzbujanja s koračno funkcijo višine \(A\) in dolžine \(B\).

import sympy as sym
t, ω_0, ω_d, d_x_0, δ, A, B, τ = sym.symbols('t, \\omega_0, \\omega_d, d_x_0, \\delta, A, B, \\tau', real=True)
π = sym.pi
podatki = {ω_0: 2, d_x_0: 1, δ: 0.2, A: 0.5, B: 10., t0: 5}

ω_d = ω_0 * sym.sqrt(1-δ**2)
ovojnica = d_x_0/ω_d * sym.exp(-δ*ω_0*t)
h = ovojnica * sym.sin(ω_d*t)*sym.Heaviside(t)

f = A*(sym.Heaviside(t)-sym.Heaviside(t-B))

fh = (f.subs(t, τ)*h.subs(t, t-τ)).subs(podatki)
fh = sym.integrate(fh, (τ, 0, t))

p1 = sym.plot(h.subs(podatki), (t,-5,20), line_color='C0', xlabel='$t$ [s]', ylabel='', title='Časovna domena',
             label='$h(t)$', show=False)
p2 = sym.plot(f.subs(podatki), (t,-5,20), line_color='C1', 
             label=f'$f(t)$, $A$={A.subs(podatki):3.1f}, $B$={B.subs(podatki):3.1f}', show=False)
p3 = sym.plot(fh, (t,-5,20), line_color='C2', 
             label='$f(t)*h(t)$', show=False)
p1.extend(p2)
p1.extend(p3)
p1.legend = True
p1.show()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 4
      2 t, ω_0, ω_d, d_x_0, δ, A, B, τ = sym.symbols('t, \\omega_0, \\omega_d, d_x_0, \\delta, A, B, \\tau', real=True)
      3 π = sym.pi
----> 4 podatki = {ω_0: 2, d_x_0: 1, δ: 0.2, A: 0.5, B: 10., t0: 5}
      6 ω_d = ω_0 * sym.sqrt(1-δ**2)
      7 ovojnica = d_x_0/ω_d * sym.exp(-δ*ω_0*t)

NameError: name 't0' is not defined

In še malo bolj podrobna slika:

import sympy as sym
t, ω_0, ω_d, d_x_0, δ, A, B, τ, t0 = sym.symbols('t, \\omega_0, \\omega_d, d_x_0, \\delta, A, B, \\tau, t0', real=True)
π = sym.pi
podatki = {ω_0: 2, d_x_0: 1, δ: 0.2, A: 0.5, B: 10., t0: 5}

ω_d = ω_0 * sym.sqrt(1-δ**2)
ovojnica = d_x_0/ω_d * sym.exp(-δ*ω_0*t)
h = ovojnica * sym.sin(ω_d*t)*sym.Heaviside(t)

f = A*(sym.Heaviside(t)-sym.Heaviside(t-B))

fh = (f.subs(t, τ)*h.subs(t, t-τ)).subs(podatki)
fh = sym.integrate(fh, (τ, -0, t))

x_array = np.linspace(0, podatki[t0], 100)
f_array = sy.lambdify(τ, f.subs(t, τ).subs(t, t0).subs(podatki)*h.subs(t, t-τ).subs(t, t0).subs(podatki))(x_array)

p1 = sym.plot(h.subs(t, t-τ).subs(t, t0).subs(podatki), (τ,-2,5), line_color='C0', 
              xlabel='$t$ [s]', ylabel='', title='Časovna domena',
             label='$h(t-\\tau)$', show=False,
             fill={'x': x_array,
                   'y1': f_array,
                   'color':'red', 'alpha':0.5})
p2 = sym.plot(f.subs(t, τ).subs(t, t0).subs(podatki), (τ,-2,5), line_color='C1', 
             label=f'$f(\\tau)$, $A$={A.subs(podatki):3.1f}, $B$={B.subs(podatki):3.1f}', show=False)
p5 = sym.plot_parametric((t0,fh.subs(t, t0)), (t0,4.95,5.05), line_color='r', 
             label=r'Vrednost $\int_{-\infty}^{t='+str(podatki[t0])+r's} f(\tau)\,h(t-\tau)\,d\tau$', show=False)
p3 = sym.plot(f.subs(t, τ).subs(t, t0).subs(podatki)*h.subs(t, t-τ).subs(t, t0).subs(podatki), (τ,-2,5), line_color='C5', 
             label=f'Podintegralska funkcija: $f(\\tau)\\,h(t-\\tau)$, $A$={A.subs(podatki):3.1f}, $B$={B.subs(podatki):3.1f}', show=False)
p4 = sym.plot(fh, (t,-1,8), line_color='C2', 
             label='$f(t)*h(t)$', show=False, nr_of_points=200)
p1.extend(p2)
p1.extend(p3)
p1.extend(p4)
p1.extend(p5)
p1.legend = True
p1.show()
../_images/8a02764871415f4441cde3b8b655f309068a278e503dd02a686ee3fc82c50568.png

Iz slike zgoraj se vidi, zakaj \(h(\tau)\) včasih poimenujemo tudi dinamični spomin sistema: \(h(t-\tau)\) namreč poskrbi, da na rešitev pri času \(t\) vpliva vse dogajanje do časa \(t\) (glejte podintegralsko funkcijo \(f(\tau)\,h(t-\tau)\) zgoraj). Iz istega razloga \(\tau\) poimenujemo spominska spremenljivka.

Opomba

Opomba: če je sistem ne-vzročen in je odziv odvisen tudi od prihodnjega vzbujanja vnosa, se meje integracije razširijo in komutativna lastnost konvolucije je:

\[ x(t) = \int_{-\infty}^{+\infty}\,f(t-\tau)\,h(\tau)\,\textrm{d} \tau = \int_{-\infty}^{+\infty}\,f(\tau)\,h(t-\tau)\,\textrm{d} \tau. \]

Povezava med frekvenčno in impulzno prenosno funkcijo#

V prejšnjih poglavjih sta bili predstavljeni frekvenčna prenosna funkcija \(\alpha(\omega)\) in impulzna prenosna funkcija \(h(t)\); prva je definirana v frekvenčni, druga v časovni domeni. V tem poglavju bomo dokazali, da \(h(t)\) in \(\alpha(\omega)\) tvorita par Fourierove transformacije.

Fourierjeva transformacija enotskega impulza \(f_1(t)=\delta(t)\):

\[ F_1(\omega)=\mathcal{F}\{f_1(t)\} =\int_{-\infty}^{\infty}f_1(t)\,\mathrm{e}^{- \mathrm{i}\,\omega\,t}\,\textrm{d} t =1 \]

in posledično enotni impulz vzbuja vse frekvence enako.

Enotski impulz prek frekvenčne prenosne funkcije \(\alpha(\omega)\) vodi v odziv \(X_1(\omega)\):

\[ X_1(\omega)=\alpha(\omega)\,\underbrace{F_1(\omega)}_{=1}=\alpha(\omega). \]

Da naredimo prehod v časovno domeno, uporabimo inverzno Fourierovo transformacijo:

\[\begin{split} \begin{split} x_1(t)&=\mathcal{F}^{-1}\{X_1(\omega)\}\\ &=\mathcal{F}^{-1}\{\alpha(\omega)\}. \end{split} \end{split}\]

Po drugi strani pa je odziv \(x_1(t)\) na enotski impulz definiran tudi z impulzno prenosno funkcijo \(x_1(t)=h(t)\), sledi:

Opomba

Impulzna prenosna funkcija in frekvenčna prenosna funkcija sta Fourierov par:

\[ h(t)=\mathcal{F}^{-1}\{\alpha(\omega)\}\qquad\textrm{ali}\qquad% \alpha(\omega)=\mathcal{F}\{h(t)\}. \]

Opomba: Fourierov par označuje dve spremenljivki, ena v časovni, druga v frekvenčni domeni, ki sta povezani prek Fourierove transformacije.