In [1]:
%%html
<style type="text/css">
input.good:checked + label {color: green}
input.bad:checked + label {color: red}
input.good:checked + label::after {color: green; content: ' Õige vastus!'}
input.bad:checked + label::after {color: red; content: ' Vale vastus!'}
</style>

Üldine $N$-mõõtmeline lineaarne süsteem

Käesolev tööleht näitab kuidas määrata üldise $N$-mõõtmelise lineaarse süsteemi püsipunkti stabiilsust, ja kuidas lahendada liikumisvõrrandeid, et leida süsteemi trajektoore.

Omaväärtused ja püsipunkti stabiilsus

Üldisel lineaarsel süsteemil $\underline{\dot{x}} = \underline{\underline{M}} \cdot \underline{x}$ on alati püsipunkt $\underline{x} = (0, \ldots, 0)$. Püsipunkti stabiilsus sõltub maatriksi $\underline{\underline{M}}$ omaväärtustest, täpsemalt nende reaalosadest. Lineaarse süsteemi püsipunkti stabiilsuse määramisel kehtivad sellised reeglid:

  • Kui kõikide omaväärtuste reaalosad on negatiivsed, siis püsipunkt on atraktor.
  • Kui kõikide omaväärtuste reaalosad on positiivsed, siis püsipunkt on repeller.
  • Kui on nii omaväärtused, mille reaalosad on positiivsed, kui omaväärtused, mille reaalosad on negatiivsed, siis püsipunkt on sadul.

Kui kõikide omaväärtuste reaalosad on nullist erinevad, siis nimetatakse püsipunkti hüperboolseks.

Üldine lahend

Lineaarse dünaamilise süsteemi lahendeid saab leida maatriksi eksponendi abil. Kui algtingimus on $\underline{x}(0) = \underline{x}_0$, siis on algtingimusele vastav lahend

$$\underline{x}(t) = e^{\underline{\underline{M}}t} \cdot \underline{x}_0$$

Maatriksi eksponent on defineeritud Taylori reana:

$$e^{\underline{\underline{M}}t} = \sum_{k = 0}^{\infty}\frac{(\underline{\underline{M}}t)^k}{k!}$$

Eksponendi saab kõige paremini arvutada baasil, milles maatriks $\underline{\underline{M}}$ on Jordani normaalkujul.

Jordani plokid ja lineaarse süsteemi normaalkuju

Igat maatriksit $\underline{\underline{M}}$ võib viia Jordani normaalkujju lineaarse teisenduse abil, see tähendab et on olemas (mitteunikaalne) pööratav maatriks $\underline{\underline{P}}$ nii, et kehtib $\underline{\underline{M}} = \underline{\underline{P}} \cdot \underline{\underline{J}} \cdot \underline{\underline{P}}^{-1}$, ja $\underline{\underline{J}}$ on kujul

$$\underline{\underline{J}} = \begin{pmatrix} \lambda_1 & 1 & & & & & & & & & & & \\ & \ddots & \ddots & & & & & & & & & & \\ & & \ddots & 1 & & & & & & & & & \\ & & & \lambda_1 & & & & & & & & & \\ & & & & \lambda_2 & 1 & & & & & & & \\ & & & & & \ddots & \ddots & & & & & & \\ & & & & & & \ddots & 1 & & & & & \\ & & & & & & & \lambda_2 & & & & & \\ & & & & & & & & \ddots & & & & & \\ & & & & & & & & & \lambda_k & 1 & & \\ & & & & & & & & & & \ddots & \ddots & \\ & & & & & & & & & & & \ddots & 1 \\ & & & & & & & & & & & & \lambda_k \\ \end{pmatrix}.$$

Siin $\lambda_1, \ldots, \lambda_k$ on maatriksi $\underline{\underline{M}}$ omaväärtused. Jordani normaalkuju koosneb plokidest

$$\underline{\underline{J}}_i = \begin{pmatrix} \lambda_i & 1 & & \\ & \ddots & \ddots & \\ & & \ddots & 1 \\ & & & \lambda_i \end{pmatrix},$$

mille suurus vastab omaväärtuse $\lambda_i$ kordsusele. Jordani maatriksi eksponent on samuti plokikujuline, ja iga plokk on vastava Jordani ploki eksponent. Selle eksponendi leidmiseks kirjutame Jordani ploki kujul $\underline{\underline{J}}_i = \underline{\underline{D}}_i + \underline{\underline{N}}_i$, kus

$$\underline{\underline{D}}_i = \begin{pmatrix} \lambda_i & & \\ & \ddots & \\ & & \lambda_i \end{pmatrix}, \quad \underline{\underline{N}}_i = \begin{pmatrix} 0 & 1 & & \\ & \ddots & \ddots & \\ & & \ddots & 1 \\ & & & 0 \end{pmatrix}.$$

Kuna need maatriksid kommuteeruvad, $\underline{\underline{D}}_i\underline{\underline{N}}_i = \underline{\underline{N}}_i\underline{\underline{D}}_i$, kehtib valem

$$e^{(\underline{\underline{D}}_i + \underline{\underline{N}}_i)t} = e^{\underline{\underline{D}}_it}e^{\underline{\underline{N}}_it}.$$

Nende maatriksi eksponendid on antud valemitega

$$e^{\underline{\underline{D}}_it} = \begin{pmatrix} e^{\lambda_it} & & \\ & \ddots & \\ & & e^{\lambda_it} \end{pmatrix}, \quad e^{\underline{\underline{N}}_it} = \begin{pmatrix} 1 & t & \cdots & \frac{t^{p_i - 1}}{(p_i - 1)!}\\ & \ddots & \ddots & \vdots \\ & & \ddots & t \\ & & & 1 \end{pmatrix},$$

kus $p_i$ on Jordani ploki $\underline{\underline{J}}_i$ suurus. Sellest järeldub, et Jordani maatriksi eksponendi kohta kehtib valem

$$e^{\underline{\underline{J}}t} = \begin{pmatrix} e^{\lambda_1t} & e^{\lambda_1t}t & \cdots & \frac{e^{\lambda_1t}t^{p_1 - 1}}{(p_1 - 1)!} & & & & & & & & & \\ & \ddots & \ddots & \vdots & & & & & & & & & \\ & & \ddots & e^{\lambda_1t}t & & & & & & & & & \\ & & & e^{\lambda_1t} & & & & & & & & & \\ & & & & e^{\lambda_2} & e^{\lambda_2}t & \cdots & \frac{e^{\lambda_2t}t^{p_2 - 1}}{(p_2 - 1)!} & & & & & \\ & & & & & \ddots & \ddots & \vdots & & & & & \\ & & & & & & \ddots & e^{\lambda_2}t & & & & & \\ & & & & & & & e^{\lambda_2} & & & & & \\ & & & & & & & & \ddots & & & & & \\ & & & & & & & & & e^{\lambda_k} & e^{\lambda_k}t & \cdots & \frac{e^{\lambda_kt}t^{p_k - 1}}{(p_k - 1)!} \\ & & & & & & & & & & \ddots & \ddots & \vdots \\ & & & & & & & & & & & \ddots & e^{\lambda_k}t \\ & & & & & & & & & & & & e^{\lambda_k} \\ \end{pmatrix}.$$

Lõpuks leiame, et maatriksi $\underline{\underline{M}}$ on

$$\begin{split} e^{\underline{\underline{M}}t} &= \sum_{j = 0}^{\infty}\frac{(\underline{\underline{M}}t)^j}{j!}\\ &= \sum_{j = 0}^{\infty}\frac{(\underline{\underline{P}} \cdot \underline{\underline{J}} \cdot \underline{\underline{P}}^{-1}t)^j}{j!}\\ &= \sum_{j = 0}^{\infty}\underline{\underline{P}} \cdot \frac{(\underline{\underline{J}}t)^j}{j!} \cdot \underline{\underline{P}}^{-1}\\ &= \underline{\underline{P}} \cdot \left(\sum_{j = 0}^{\infty}\frac{(\underline{\underline{J}}t)^j}{j!}\right) \cdot \underline{\underline{P}}^{-1}\\ &=\underline{\underline{P}} \cdot e^{\underline{\underline{J}}t} \cdot \underline{\underline{P}}^{-1}. \end{split}$$

Interaktiivne näide: maatriksi eksponent

Järgmisena rakendame oma teadmisi praktikale. Alljärgnev kood näitab, kuidas saab maatriksi eksponendi leida SymPy abil, ja sellega ka lineaarse süsteemi trajektoore leida. Kõigepealt on meil selleks SymPy paketti vaja.

In [2]:
import sympy as sp
from sympy.interactive.printing import init_printing
init_printing(use_unicode=True)

Siis defineerime sümboleid: aega tähistame $t$-ga, ja algtingimus on $(x_{10}, x_{20}, x_{30})$.

In [3]:
t = sp.symbols('t', real=True)
x0 = sp.Matrix(3, 1, sp.symbols('x_10 x_20 x_30', real=True))

Näide #1

Esimene näide on süsteem, mida kirjeldab dünaamika

$$\begin{pmatrix} \dot{x}_1\\ \dot{x}_2\\ \dot{x}_3 \end{pmatrix} = \begin{pmatrix} 1 & -\frac{2}{3} & 0\\ 2 & 3 & 1\\ 0 & -\frac{2}{3} & 1\\ \end{pmatrix} \cdot \begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix}.$$

Maatriks $\underline{\underline{M}}$ on seega selline:

In [4]:
M = sp.Matrix([[1, -sp.Rational(2, 3), 0], [2, 3, 1], [0, -sp.Rational(2, 3), 1]])
M
Out[4]:
$\displaystyle \left[\begin{matrix}1 & - \frac{2}{3} & 0\\2 & 3 & 1\\0 & - \frac{2}{3} & 1\end{matrix}\right]$

Vaatame, millised on maatriksi omaväärtused.

In [5]:
ev = M.eigenvals()
ev
Out[5]:
$\displaystyle \left\{ 1 : 1, \ 2 - i : 1, \ 2 + i : 1\right\}$

Maatriksil on üks reaalne omaväärtus $\lambda_1 = 1$ ja kaaskompleksne paar $\lambda_{2,3} = 2 \pm i$. Reaalosad on positiivsed, seega on püsipunkt $(0, 0, 0)$ repeller. Iga omaväärtuse kordsus on 1, seega kaasneb iga omaväärtusega üks omavektor.

In [6]:
[(M - y * sp.eye(3)).nullspace()[0].applyfunc(sp.nsimplify) for y in ev]
Out[6]:
$\displaystyle \left[ \left[\begin{matrix}- \frac{1}{2}\\0\\1\end{matrix}\right], \ \left[\begin{matrix}1\\- \frac{3}{2} + \frac{3 i}{2}\\1\end{matrix}\right], \ \left[\begin{matrix}1\\- \frac{3}{2} - \frac{3 i}{2}\\1\end{matrix}\right]\right]$

Omavektorid moodustavad teisendusmaatriksi $\underline{\underline{P}}$.

In [7]:
P, jc = M.jordan_cells()
P = P.applyfunc(sp.nsimplify)
P
Out[7]:
$\displaystyle \left[\begin{matrix}- \frac{1}{2} & 1 & 1\\0 & - \frac{3}{2} + \frac{3 i}{2} & - \frac{3}{2} - \frac{3 i}{2}\\1 & 1 & 1\end{matrix}\right]$

Selle abil leiame ka Jordani normaalkuju $\underline{\underline{J}}$.

In [8]:
J = (P.inv() * M * P).applyfunc(sp.nsimplify)
J
Out[8]:
$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 2 - i & 0\\0 & 0 & 2 + i\end{matrix}\right]$

Nüüd on ka lihtne leida selle maatriksi eksponenti $e^{\underline{\underline{J}}t}$.

In [9]:
eJ = (t * J).exp()
eJ
Out[9]:
$\displaystyle \left[\begin{matrix}e^{t} & 0 & 0\\0 & e^{t \left(2 - i\right)} & 0\\0 & 0 & e^{t \left(2 + i\right)}\end{matrix}\right]$

Lõpuks kasutame lineaarset teisendust, et leida maatriksi eksponenti $e^{\underline{\underline{M}}t}$.

In [10]:
eM = (P * eJ * P.inv()).applyfunc(lambda y: sp.expand(sp.expand_trig(sp.simplify(sp.simplify(y).as_real_imag()[0]))))
eM
Out[10]:
$\displaystyle \left[\begin{matrix}- \frac{2 e^{2 t} \sin{\left(t \right)}}{3} + \frac{2 e^{2 t} \cos{\left(t \right)}}{3} + \frac{e^{t}}{3} & - \frac{2 e^{2 t} \sin{\left(t \right)}}{3} & - \frac{e^{2 t} \sin{\left(t \right)}}{3} + \frac{e^{2 t} \cos{\left(t \right)}}{3} - \frac{e^{t}}{3}\\2 e^{2 t} \sin{\left(t \right)} & e^{2 t} \sin{\left(t \right)} + e^{2 t} \cos{\left(t \right)} & e^{2 t} \sin{\left(t \right)}\\- \frac{2 e^{2 t} \sin{\left(t \right)}}{3} + \frac{2 e^{2 t} \cos{\left(t \right)}}{3} - \frac{2 e^{t}}{3} & - \frac{2 e^{2 t} \sin{\left(t \right)}}{3} & - \frac{e^{2 t} \sin{\left(t \right)}}{3} + \frac{e^{2 t} \cos{\left(t \right)}}{3} + \frac{2 e^{t}}{3}\end{matrix}\right]$

Algtingimusele $(x_{10}, x_{20}, x_{30})$ vastav lahend on siis selline.

In [11]:
x = (eM * x0).applyfunc(sp.simplify)
x
Out[11]:
$\displaystyle \left[\begin{matrix}\frac{\left(2 \sqrt{2} x_{10} e^{t} \cos{\left(t + \frac{\pi}{4} \right)} + x_{10} - 2 x_{20} e^{t} \sin{\left(t \right)} + \sqrt{2} x_{30} e^{t} \cos{\left(t + \frac{\pi}{4} \right)} - x_{30}\right) e^{t}}{3}\\\left(2 x_{10} \sin{\left(t \right)} + \sqrt{2} x_{20} \sin{\left(t + \frac{\pi}{4} \right)} + x_{30} \sin{\left(t \right)}\right) e^{2 t}\\\frac{\left(- 2 x_{10} \left(- \sqrt{2} e^{t} \cos{\left(t + \frac{\pi}{4} \right)} + 1\right) - 2 x_{20} e^{t} \sin{\left(t \right)} + x_{30} \left(\sqrt{2} e^{t} \cos{\left(t + \frac{\pi}{4} \right)} + 2\right)\right) e^{t}}{3}\end{matrix}\right]$

Igaks juhuks kontrollime, kas dünaamiline võrrand an rahuldatud.

In [12]:
(x.applyfunc(lambda y: sp.diff(y, t)) - M * x).applyfunc(sp.simplify)
Out[12]:
$\displaystyle \left[\begin{matrix}0\\0\\0\end{matrix}\right]$

Näide #2

Teine näide on süsteem, mida kirjeldab dünaamika

$$\begin{pmatrix} \dot{x}_1\\ \dot{x}_2\\ \dot{x}_3 \end{pmatrix} = \begin{pmatrix} 1 & 2 & 0\\ 1 & 1 & 2\\ 0 & -1 & 1\\ \end{pmatrix} \cdot \begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix}.$$

Maatriks $\underline{\underline{M}}$ on seega selline:

In [13]:
M = sp.Matrix([[1, 2, 0], [1, 1, 2], [0, -1, 1]])
M
Out[13]:
$\displaystyle \left[\begin{matrix}1 & 2 & 0\\1 & 1 & 2\\0 & -1 & 1\end{matrix}\right]$

Vaatame, millised on maatriksi omaväärtused.

In [14]:
ev = M.eigenvals()
ev
Out[14]:
$\displaystyle \left\{ 1 : 3\right\}$

Maatriksil on üks reaalne omaväärtus $\lambda_1 = 1$, mille kordsus on 3. Reaalosa on positiivne, seega on püsipunkt $(0, 0, 0)$ repeller. Sellega kaasneb üks omavektor.

In [15]:
[(M - y * sp.eye(3)).nullspace()[0].applyfunc(sp.nsimplify) for y in ev]
Out[15]:
$\displaystyle \left[ \left[\begin{matrix}-2\\0\\1\end{matrix}\right]\right]$

Omavektor ja üldistatud omavektorid moodustavad teisendusmaatriksi $\underline{\underline{P}}$.

In [16]:
P, jc = M.jordan_cells()
P = P.applyfunc(sp.nsimplify)
P
Out[16]:
$\displaystyle \left[\begin{matrix}2 & 0 & 1\\0 & 1 & 0\\-1 & 0 & 0\end{matrix}\right]$

Selle abil leiame ka Jordani normaalkuju $\underline{\underline{J}}$.

In [17]:
J = (P.inv() * M * P).applyfunc(sp.nsimplify)
J
Out[17]:
$\displaystyle \left[\begin{matrix}1 & 1 & 0\\0 & 1 & 1\\0 & 0 & 1\end{matrix}\right]$

Nüüd on ka lihtne leida selle maatriksi eksponenti $e^{\underline{\underline{J}}t}$.

In [18]:
eJ = (t * J).exp()
eJ
Out[18]:
$\displaystyle \left[\begin{matrix}e^{t} & t e^{t} & \frac{t^{2} e^{t}}{2}\\0 & e^{t} & t e^{t}\\0 & 0 & e^{t}\end{matrix}\right]$

Lõpuks kasutame lineaarset teisendust, et leida maatriksi eksponenti $e^{\underline{\underline{M}}t}$.

In [19]:
eM = (P * eJ * P.inv()).applyfunc(sp.simplify)
eM
Out[19]:
$\displaystyle \left[\begin{matrix}\left(t^{2} + 1\right) e^{t} & 2 t e^{t} & 2 t^{2} e^{t}\\t e^{t} & e^{t} & 2 t e^{t}\\- \frac{t^{2} e^{t}}{2} & - t e^{t} & \left(1 - t^{2}\right) e^{t}\end{matrix}\right]$

Algtingimusele $(x_{10}, x_{20}, x_{30})$ vastav lahend on siis selline.

In [20]:
x = eM * x0
x
Out[20]:
$\displaystyle \left[\begin{matrix}2 t^{2} x_{30} e^{t} + 2 t x_{20} e^{t} + x_{10} \left(t^{2} + 1\right) e^{t}\\t x_{10} e^{t} + 2 t x_{30} e^{t} + x_{20} e^{t}\\- \frac{t^{2} x_{10} e^{t}}{2} - t x_{20} e^{t} + x_{30} \left(1 - t^{2}\right) e^{t}\end{matrix}\right]$

Igaks juhuks kontrollime, kas dünaamiline võrrand an rahuldatud.

In [21]:
(x.applyfunc(lambda y: sp.diff(y, t)) - M * x).applyfunc(sp.simplify)
Out[21]:
$\displaystyle \left[\begin{matrix}0\\0\\0\end{matrix}\right]$

Küsimused

Hüperboolne püsipunkt

Millised on Jacobi maatriksi omaväärtuste reaalosad hüperboolse püsipunkti juures?

Polünoomid maatriksi eksponendis

Millisel juhul sisaldab maatriksi eksponent $e^{\underline{\underline{M}}t}$ polünoome ($t, t^2, \ldots$)?

Püsipunkti tüüp

Milline on järgneva süsteemi püsipunkt $\underline{x}_* = (0, 0, 0)$?

$$\begin{pmatrix} \dot{x}_1\\ \dot{x}_2\\ \dot{x}_3 \end{pmatrix} = \begin{pmatrix} -5 & -4 & 9 \\ 3 & 2 & -3 \\ -4 & -4 & 8 \end{pmatrix} \cdot \begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix}$$