SharedRovniceNB.sagewsOpen in CoCalc

Náhodné procesy:

riešenie sústavy diferenciálnych rovníc

metódou vlastných čísel

\displaystyle \frac{dp_{0}(t)}{d t} = -2 \lambda p_{0}(t) + \mu p_{1}(t) \\ \displaystyle\frac{dp_{1}(t)}{d t} = 2 \lambda p_{0}(t) - {(\lambda + \mu)} p_{1}(t) + 2 \mu p_{2}(t) \\ \displaystyle\frac{dp_{2}(t)}{d t} = \lambda p_{1}(t) - 2 \mu p_{2}(t)

p(0)=1,\, p_1(0)=p_2(0)=0, \, \lambda, \mu \in \mathbb{R}^+

V teórii diferenciálnych rovníc možno tento systém riešiť:

  • pomocou Laplaceovej transformácie (viacrozmerná analógia s 1D LT)
  • pomocou vlastných hodnôt (viacrozmerná analógia s riešením rovnice tretieho rádu)
  • pomocou exponenciálnej matice (viacrozmerná analógia s riešením rovnice prvého rádu)

Ref.
JOYNER, D. a HAMPTON, M., 2012. Introduction to Differential Equations Using Sage. Baltimore: Johns Hopkins University Press. ISBN 978-1-4214-0637-4.
ABELL, M.L.L. a BRASELTON, J.P., 2016. Differential Equations with Mathematica, Fourth Edition. 4 edition. S.l.: Academic Press. ISBN 978-0-12-804776-7.
KLUVÁNEK, I., MIŠÍK, L. a ŠVEC, M., 1961. Matematika II. 3. Bratislava: Alfa.

Riešenie pomocou Sage

#Zadanie parametrov lambda, mu
%var l,m
assume(l>0,m>0)
#Zadanie hľadaných funkcií
p0(t)=function('p0')(t)
p1(t)=function('p1')(t)
p2(t)=function('p2')(t)
funkcie = [p0,p1,p2]
#Zadanie sústavy diferenciálnych rovníc
dr0 = diff(p0,t) == -2*l*p0+m*p1
dr1 = diff(p1,t) == -(l+m)*p1+2*l*p0+2*m*p2
dr2 = diff(p2,t) == l*p1-2*m*p2
system = [dr0,dr1,dr2]
#množina riešení sústavy
mries = desolve_system(system,funkcie,ivar=t, ics=[0,1,0,0])
#výpis riešení
for riesenie in mries:
    show(riesenie)
\displaystyle p_{0}\left(t\right) = \frac{2 \, l m e^{\left(-{\left(l + m\right)} t\right)}}{l^{2} + 2 \, l m + m^{2}} + \frac{l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{l^{2} + 2 \, l m + m^{2}} + \frac{m^{2}}{l^{2} + 2 \, l m + m^{2}}
\displaystyle p_{1}\left(t\right) = -\frac{2 \, l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{l^{2} + 2 \, l m + m^{2}} + \frac{2 \, l m}{l^{2} + 2 \, l m + m^{2}} + \frac{2 \, {\left(l^{2} - l m\right)} e^{\left(-{\left(l + m\right)} t\right)}}{l^{2} + 2 \, l m + m^{2}}
\displaystyle p_{2}\left(t\right) = -\frac{2 \, l^{2} e^{\left(-{\left(l + m\right)} t\right)}}{l^{2} + 2 \, l m + m^{2}} + \frac{l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{l^{2} + 2 \, l m + m^{2}} + \frac{l^{2}}{l^{2} + 2 \, l m + m^{2}}
#výpis riešení v zjednodušenom tvare
for riesenie in mries:
    show(riesenie.factor())
\displaystyle p_{0}\left(t\right) = \frac{{\left(m e^{\left(l t + m t\right)} + l\right)}^{2} e^{\left(-2 \, l t - 2 \, m t\right)}}{{\left(l + m\right)}^{2}}
\displaystyle p_{1}\left(t\right) = \frac{2 \, {\left(m e^{\left(l t + m t\right)} + l\right)} l {\left(e^{\left(l t + m t\right)} - 1\right)} e^{\left(-2 \, l t - 2 \, m t\right)}}{{\left(l + m\right)}^{2}}
\displaystyle p_{2}\left(t\right) = \frac{l^{2} {\left(e^{\left(l t + m t\right)} - 1\right)}^{2} e^{\left(-2 \, l t - 2 \, m t\right)}}{{\left(l + m\right)}^{2}}
#zápis rovníc v texu
text =''
for rovnica in system:
    text+=latex(rovnica)
text
t \ {\mapsto}\ \frac{\partial}{\partial t}p_{0}\left(t\right) = -2 \, l p_{0}\left(t\right) + m p_{1}\left(t\right) t \ {\mapsto}\ \frac{\partial}{\partial t}p_{1}\left(t\right) = 2 \, l p_{0}\left(t\right) - {\left(l + m\right)} p_{1}\left(t\right) + 2 \, m p_{2}\left(t\right) t \ {\mapsto}\ \frac{\partial}{\partial t}p_{2}\left(t\right) = l p_{1}\left(t\right) - 2 \, m p_{2}\left(t\right)

Riešenie krok po kroku

Metóda vlastných čísel a charakteristického polynómu
analogická metóde riešenia lineárnej diferenciálnej rovnice
s konštantnými koeficientami

1. krok - zápis do maticového tvaru

Systém diferenciálnych rovníc prepíšeme do maticového tvaru

\left(\begin{array}{r} p_{0}'\left(t\right) \\ p_{1}'\left(t\right) \\ p_{2}'\left(t\right) \end{array}\right) =\left(\begin{array}{rrr} -2 \, \lambda & \mu & 0 \\ 2 \, \lambda & -\lambda - \mu & 2 \, \mu \\ 0 & \lambda & -2 \, \mu \end{array}\right)\left(\begin{array}{r} p_{0}\left(t\right) \\ p_{1}\left(t\right) \\ p_{2}\left(t\right) \end{array}\right)

\phantom{xx}\vec{p}'(t)\hspace{0.5cm}= \hspace{2cm}\mathrm{A}\hspace{3cm}\vec{p}(t)

2. krok - hľadáme fundamentálne riešenie s exponenciálou

Riešenie hľadáme v tvare:
\vec{p}(t) = \vec{v}e^{rt}, \, r\in \mathbb{C}, \vec{v} \in \mathbb{C}^3

Dosadením do sústavy musí platiť pre \vec{v},r :
\vec{v}re^{rt}=\mathrm{A}\vec{v}e^{rt}
\hspace{0.5cm}r\vec{v}=\mathrm{A}\vec{v}

Skalár r je vlastným číslom \mathrm{A} a vektor \vec{v} nenulovým vlastným vektorom \mathrm{A} .

3. krok - nájdeme všetky fundamentálne riešenia

t.j. pomocou vlastných čísel a vektorov \mathrm{A}

Vlastné čísla a nenulové vektory sú netriviálnym riešením:
\hspace{0.5cm}r\vec{v}=\mathrm{A}\vec{v}

Prepísaním tejto sústavy do tvaru:
\hspace{0.5cm}\mathrm{A}\vec{v}-r\mathrm{I}\vec{v}=\vec{0}
\hspace{0.3cm}\left(\mathrm{A}-r\mathrm{I}\right)\vec{v}=\vec{0}

Netriviálne riešenia dostaneme len v prípade singulárnej matice \left(\mathrm{A}-r\mathrm{I}\right) ,
t.j. jej determinant musí byť nula.

#Zadanie matice A
A = matrix(SR, 3, 3, [-2*l,m,0,2*l,-l-m,2*m,0,l,-2*m])
show(A)
\displaystyle \left(\begin{array}{rrr} -2 \, l & m & 0 \\ 2 \, l & -l - m & 2 \, m \\ 0 & l & -2 \, m \end{array}\right)
#vytvorenie matice A-rI
%var r
II = identity_matrix(3)
show(A-r*II)
\displaystyle \left(\begin{array}{rrr} -2 \, l - r & m & 0 \\ 2 \, l & -l - m - r & 2 \, m \\ 0 & l & -2 \, m - r \end{array}\right)
# determinant A-rI
poly=(A-r*II).det().factor()
poly
-(2*l + 2*m + r)*(l + m + r)*r
#vlastné čísla ako riešenie rovnice det(A-rI) = 0
solve(poly, r)
[r == -2*l - 2*m, r == -l - m, r == 0]
#priamy výpočet vlastných hodnôt v Sagi
vhodnoty = A.eigenvalues()
vhodnoty
[-2*l - 2*m, -l - m, 0]
#vlastný vektor k r=0, riešenie (A-0I)v=0
# Gaussova metóda - trojuholníkový tvar
(A-0*II).rref()
[ 1 0 -m^2/l^2] [ 0 1 -2*m/l] [ 0 0 0]
#vlastný vektor
r = [0]
v = [vector([m^2/l^2,2*m/l,1])]
show(v[0])
(A-r[0]*II)*v[0]
_.canonicalize_radical()
\displaystyle \left(\frac{m^{2}}{l^{2}},\,\frac{2 \, m}{l},\,1\right)
(0, -2*(l + m)*m/l + 2*m + 2*m^2/l, 0) (0, 0, 0)
#vlastný vektor k r=-l-m
(A+(l+m)*II).rref()
[ 1 0 m/l] [ 0 1 (l - m)/l] [ 0 0 0]
#vlastný vektor
r += [-(l+m)]
v += [vector([-m/l,(m-l)/l,1])]
show(v[1])
(A-r[1]*II)*v[1]
\displaystyle \left(-\frac{m}{l},\,-\frac{l - m}{l},\,1\right)
(0, 0, 0)
#vlastný vektor k r=-2*(l+m)
(A+2*(l+m)*II).rref()
[ 1 0 -1] [ 0 1 2] [ 0 0 0]
#vlastný vektor
r +=[-2*(l+m)]
v +=[vector([1,-2,1])]
show(v[2])
(A-r[2]*II)*v[2]
\displaystyle \left(1,\,-2,\,1\right)
(0, 0, 0)
# priamy výpočet vlastných vektorov v Sagi - ten vypočítava druhú množinu vektorov
vlastneA = A.eigenmatrix_right()
vecA = vlastneA[1]
show(vecA)
\displaystyle \left(\begin{array}{rrr} 1 & 1 & 1 \\ -2 & \frac{l - m}{m} & \frac{2 \, l}{m} \\ 1 & -\frac{l}{m} & \frac{l^{2}}{m^{2}} \end{array}\right)

5. krok - zapíšeme všeobecné riešenie pomocou fundamentálnych

\vec{p}(t) = C_1 \vec{v}_1 e^{r_1t}+C_2 \vec{v}_2 e^{r_2t}+C_3 \vec{v}_3 e^{r_3t}, C_1, C_2, C_3 \in \mathbb{R}^3

6.krok - nájdeme riešenie s počiatočnými podmienkami

t.j. nájdeme konštanty C_1, C_2, C_3

Tieto konštanty sú riešením sústavy: \left(\begin{array}{r} p_{0}\left(0\right) \\ p_{1}\left(0\right) \\ p_{2}\left(0\right) \end{array}\right) = \mathrm{V} \left(\begin{array}{r} C_1 \\ C_2 \\ C_3 \end{array}\right)

\mathrm{V}=\left(\vec{v}_1 \vec{v}_2 \vec{v}_3 \right)

# konštrukcia matice V
V = matrix(3,3,[v[0],v[1],v[2]]).T
show(V)
InvV = (V^(-1)).factor()
show(InvV)
\displaystyle \left(\begin{array}{rrr} \frac{m^{2}}{l^{2}} & -\frac{m}{l} & 1 \\ \frac{2 \, m}{l} & -\frac{l - m}{l} & -2 \\ 1 & 1 & 1 \end{array}\right)
\displaystyle \left(\begin{array}{rrr} \frac{l^{2}}{{\left(l + m\right)}^{2}} & \frac{l^{2}}{{\left(l + m\right)}^{2}} & \frac{l^{2}}{{\left(l + m\right)}^{2}} \\ -\frac{2 \, l^{2}}{{\left(l + m\right)}^{2}} & -\frac{{\left(l - m\right)} l}{{\left(l + m\right)}^{2}} & \frac{2 \, l m}{{\left(l + m\right)}^{2}} \\ \frac{l^{2}}{{\left(l + m\right)}^{2}} & -\frac{l m}{{\left(l + m\right)}^{2}} & \frac{m^{2}}{{\left(l + m\right)}^{2}} \end{array}\right)
# Výpočet konštánt C
C = InvV*vector([1,0,0])
show(C)
\displaystyle \left(\frac{l^{2}}{{\left(l + m\right)}^{2}},\,-\frac{2 \, l^{2}}{{\left(l + m\right)}^{2}},\,\frac{l^{2}}{{\left(l + m\right)}^{2}}\right)
# konštanty pre počiatočnú podmienku v Sagi
Cs = (vecA^(-1)*vector([1,0,0])).canonicalize_radical()
show(Cs)
\displaystyle \left(\frac{l^{2}}{l^{2} + 2 \, l m + m^{2}},\,\frac{2 \, l m}{l^{2} + 2 \, l m + m^{2}},\,\frac{m^{2}}{l^{2} + 2 \, l m + m^{2}}\right)
# riešenie s počiatočnou podmienkou
%var t
pt = sum(C[n]*exp(r[n]*t)*v[n] for n in [0,1,2])
for n in [0,1,2]:
    show(funkcie[n](t)==pt[n])
\displaystyle p_{0}\left(t\right) = \frac{2 \, l m e^{\left(-{\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{m^{2}}{{\left(l + m\right)}^{2}}
\displaystyle p_{1}\left(t\right) = \frac{2 \, {\left(l - m\right)} l e^{\left(-{\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} - \frac{2 \, l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{2 \, l m}{{\left(l + m\right)}^{2}}
\displaystyle p_{2}\left(t\right) = -\frac{2 \, l^{2} e^{\left(-{\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{l^{2}}{{\left(l + m\right)}^{2}}
for n in [0,1,2]:
    latex(funkcie[n](t)==pt[n])
p_{0}\left(t\right) = \frac{2 \, l m e^{\left(-{\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{m^{2}}{{\left(l + m\right)}^{2}} p_{1}\left(t\right) = \frac{2 \, {\left(l - m\right)} l e^{\left(-{\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} - \frac{2 \, l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{2 \, l m}{{\left(l + m\right)}^{2}} p_{2}\left(t\right) = -\frac{2 \, l^{2} e^{\left(-{\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{l^{2} e^{\left(-2 \, {\left(l + m\right)} t\right)}}{{\left(l + m\right)}^{2}} + \frac{l^{2}}{{\left(l + m\right)}^{2}}
#porovnanie s riešením SAGE
for riesenie in mries:
    show(riesenie.factor())
\displaystyle p_{0}\left(t\right) = \frac{{\left(m e^{\left(l t + m t\right)} + l\right)}^{2} e^{\left(-2 \, l t - 2 \, m t\right)}}{{\left(l + m\right)}^{2}}
\displaystyle p_{1}\left(t\right) = \frac{2 \, {\left(m e^{\left(l t + m t\right)} + l\right)} l {\left(e^{\left(l t + m t\right)} - 1\right)} e^{\left(-2 \, l t - 2 \, m t\right)}}{{\left(l + m\right)}^{2}}
\displaystyle p_{2}\left(t\right) = \frac{l^{2} {\left(e^{\left(l t + m t\right)} - 1\right)}^{2} e^{\left(-2 \, l t - 2 \, m t\right)}}{{\left(l + m\right)}^{2}}

Výpočet na papieri - zhrnutie