Линейные многомерные модели

Уравнение стохастического осциллятора <<	Оглавление	>> Многомерие помогает одномерию

---

$\textstyle \bullet$ Найдём решение линейных стохастических уравнений (по $\textstyle j$ — сумма):

$dx_{i}=A_{ij}\cdot (x_{j}-c_{j})\,dt+B_{ij}\cdot \delta {W_{j}}.$

Постоянный вектор $\textstyle c_{j}$ можно убрать сдвигом $\textstyle x_{j}\to x_{j}+c_{j}$. В решении делается обратный сдвиг. Поэтому будем изучать однородное уравнение, которое запишем в матричной форме:

$d\mathbf {x} =\mathbf {A} \cdot \mathbf {x} \;dt+\mathbf {B} \cdot \delta \mathbf {W} ,$

где $\textstyle \mathbf {A}$ и $\textstyle \mathbf {B}$ — не зависящие от $\textstyle \mathbf {x}$ и времени матрицы.

Для определения среднего проще всего сразу воспользоваться соотношением (6.16):

${\dot {\overline {\mathbf {x} }}}=\mathbf {A} \cdot {\overline {\mathbf {x} }}\;\;\;\;\;\;=>\;\;\;\;\;\;\;\;{\overline {\mathbf {x} }}=e^{\mathbf {A} t}\cdot \mathbf {x} _{0},$

(6.20)

где $\textstyle \mathbf {x} _{0}$ — вектор начального значения. Если мы хотим "вернуть" $\textstyle \mathbf {c}$, то потребуются две замены: $\textstyle {\overline {\mathbf {x} }}\to {\overline {\mathbf {x} }}-\mathbf {c}$ и $\textstyle \mathbf {x} _{0}\to \mathbf {x} _{0}-\mathbf {c}$.

$\textstyle \bullet$ Монотонная зависимость от $\textstyle t$ в матричной записи решения (6.20) обманчива. Рассмотрим стохастический осциллятор из предыдущего раздела:

${\begin{pmatrix}dx\\dy\\\end{pmatrix}}={\begin{pmatrix}-\lambda &-\omega \\\omega &-\lambda \\\end{pmatrix}}\cdot {\begin{pmatrix}x\\y\\\end{pmatrix}}+\sigma \cdot {\begin{pmatrix}\delta W_{x}\\\delta W_{y}\\\end{pmatrix}}.$

(6.21)

В этом случае матрицу $\textstyle \mathbf {A}$ можно разбить на сумму двух матриц:

Невозможно разобрать выражение (MathML с переходом в SVG или PNG (рекомендуется для современных браузеров и инструментов повышения доступности): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle \mathbf{A } = \begin{pmatrix} -\lambda & -\omega \\ \omega & -\lambda \\ \end{pmatrix} = \omega \cdot \mathbf{q} -\lambda \cdot \mathbf{1},\;\;\;\;\;\;где\;\;\;\; \mathbf{1}= \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}, \;\;\;\;\;\; \mathbf{q}= \begin{pmatrix} 0 & -1 \\ 1 & 0 \\ \end{pmatrix}.}

Несложно проверить, что:

$\mathbf {q} ^{2}=-\mathbf {1} ,\;\;\;\mathbf {q} ^{3}=-\mathbf {q} ,\;\;\;\mathbf {q} ^{4}=\mathbf {1} ,\;\;\;\;\mathbf {q} ^{5}=\mathbf {q} ,...$

Так как матрицы $\textstyle \mathbf {1}$ и $\textstyle \mathbf {q}$ коммутируют друг с другом ($\textstyle \mathbf {q} \cdot \mathbf {1} =\mathbf {1} \cdot \mathbf {q}$), экспонента суммы разбивается на произведение $\textstyle e^{\mathbf {A} t}=e^{-\mathbf {1} \lambda t}\cdot e^{\mathbf {q} \,\omega t}$. Раскладывая второй множитель по степеням $\textstyle t$ и учитывая аналогичное разложение для синуса и косинуса, решение можно представить в следующем виде:

${\begin{pmatrix}{\bar {x}}\\{\bar {y}}\\\end{pmatrix}}=e^{-\lambda t}\cdot {\begin{pmatrix}\cos \omega t&-\sin \omega t\\\sin \omega t&\cos \omega t\\\end{pmatrix}}{\begin{pmatrix}x_{0}\\y_{0}\\\end{pmatrix}}.$

(6.22)

Оно же выше было получено другим методом. Таким образом, "монотонная" зависимость от времени в матричных соотношениях вполне может превратиться в периодическую функцию.

$\textstyle \bullet$ Найдём более практичное, чем (6.20), представление для решения линейного уравнения. Будем его искать в виде:

${\overline {\mathbf {x} }}(t)=\mathbf {u} \,e^{at}\;\;\;\;\;\;\;\;\;\;=>\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathbf {A} \cdot \mathbf {u} =a\,\mathbf {u} .$

(6.23)

Постоянный вектор $\textstyle \mathbf {u}$ является собственным вектором матрицы $\textstyle \mathbf {A}$, а параметр "$\textstyle a$"— её собственным значением. Перенося $\textstyle (a\,\mathbf {u} )$ в левую часть, получаем систему однородных уравнений относительно $\textstyle \mathbf {u}$, которая имеет ненулевое решение, только если её детерминант равен нулю:

$\det(\mathbf {A} -a\,\mathbf {1} )=0.$

Это уравнение называется характеристическим и является полиномом $\textstyle n$-той степени по $\textstyle a$. Обычно оно имеет $\textstyle n$ различных решений $\textstyle a_{1},...,a_{n}$. Часть из них может оказаться комплексными. Для каждого из них мы решаем уравнение (6.23) и находим собственные вектора $\textstyle \mathbf {u} ^{(k)}$. Внимание! Верхний индекс — это номер собственного вектора, а не его компонента.

Теперь общее решение для среднего значения вектора переменных состояния можно записать в следующем виде:

${\overline {\mathbf {x} }}(t)=\sum _{k}\mu _{k}\,\mathbf {u} ^{(k)}\,e^{a_{k}t},\;\;\;\;\;\;\;\;\;\;\;\;\mathbf {x} _{0}=\sum _{k}\mu _{k}\,\mathbf {u} ^{(k)},$

(6.24)

где $\textstyle \mu _{k}$ — произвольные константы, выражающиеся через начальные условия $\textstyle \mathbf {x} _{0}=\mathbf {x} (0)$. Прямой подстановкой в исходное уравнение можно проверить справедливость этого решения. Действительная часть собственных значений $\textstyle a_{k}$ будет приводить к экспоненциально уменьшающимся ($\textstyle \mathbf {Re} \,a_{k}<0$) или увеличивающимся ($\textstyle \mathbf {Re} \,a_{k}>0$) решениям. Мнимая часть соответствует колебательным режимам.

Если матрица $\textstyle \mathbf {A}$ симметрична, то собственные вектора можно выбрать ортогональными: $\textstyle \mathbf {u} ^{(\alpha )}\cdot \mathbf {u} ^{*(\beta )}=\delta _{\alpha \beta }$ (звёздочка — комплексное сопряжение). В этом случае $\textstyle \mu _{k}=\mathbf {x} _{0}\cdot \mathbf {u} ^{*(k)}$.

Когда $\textstyle \mu _{k}$ выражены через $\textstyle \mathbf {x} _{0}$, можно найти явное представление экспоненты от матрицы. Действительно, из (6.20), взяв производную по компонентам начального условия, имеем $\textstyle \left[e^{\mathbf {A} t}\right]_{\alpha \beta }=\partial {\overline {x}}_{\alpha }/x_{0\beta }$. В частности, если собственные вектора ортогональны ($\textstyle \mu _{k}=\mathbf {x} _{0}\cdot \mathbf {u} ^{*(k)}$), то:

$\left[e^{\mathbf {A} t}\right]_{\alpha \beta }=\sum _{k}u_{\alpha }^{(k)}\,u_{\beta }^{*(k)}\,e^{a_{k}t}.$

(6.25)

В качестве упражнения ($\textstyle \lessdot$ H) предлагается найти $\textstyle e^{\mathbf {A} t}$ для матрицы 2x2, у которой $\textstyle A_{12}=A_{22}=0$. Необходимо это сделать прямым разложением экспоненты при помощи собственных значений.

$\textstyle \bullet$ Выразим теперь решение стохастической линейной системы через гауссовы переменные. Введём новый вектор $\textstyle \mathbf {y}$, удовлетворяющий, по лемме Ито (6.13), следующему уравнению:

$\mathbf {y} =e^{-\mathbf {A} t}\cdot \mathbf {x} \;\;\;\;\;=>\;\;\;\;\;\;d\mathbf {y} =e^{-\mathbf {A} t}\,\mathbf {B} \,\delta \mathbf {W} =\mathbf {G} (t)\,\delta \mathbf {W} .$

Матрица $\textstyle \mathbf {G} (t)=e^{-\mathbf {A} t}\,\mathbf {B}$ зависит только от времени, поэтому решение этого уравнения легко найти при помощи итерационного метода:

$y_{\mu }(t)=y_{\mu }(t_{0})+\sum _{k}G_{\mu \alpha }(t_{k})\varepsilon _{\alpha }(t_{k}){\sqrt {\Delta t}}=y_{\mu }(t_{0})+g_{\mu \alpha }\varepsilon _{\alpha }.$

Сумма независимых гауссовых чисел $\textstyle \varepsilon _{\alpha }(t_{k})$ снова пропорциональна гауссовому числу, которое удобно представить в виде суммы независимых величин $\textstyle \varepsilon _{\alpha }$ (второе равенство). Найдём значения $\textstyle g_{\mu \alpha }$. Для этого вычислим среднее от $\textstyle \left\langle {\bigl (}y(t)-y(t_{0}))_{\mu }(y(t)-y(t_{0}){\bigr )}_{\nu }\right\rangle$:

$g_{\mu \alpha }g_{\nu \beta }\left\langle \varepsilon _{\alpha }\varepsilon _{\beta }\right\rangle =\sum _{k,l}G_{\mu \alpha }(t_{k})G_{\nu \beta }(t_{l})\left\langle \varepsilon _{\alpha }(t_{k})\varepsilon _{\beta }(t_{l})\right\rangle \Delta t.$

Учитывая независимость случайных величин $\textstyle \left\langle \varepsilon _{\alpha }(t_{k})\varepsilon _{\beta }(t_{l})\right\rangle =\delta _{\alpha ,\beta }\delta _{k,l}$ и $\textstyle \left\langle \varepsilon _{\alpha }\varepsilon _{\beta }\right\rangle =\delta _{\alpha ,\beta }$, а также переходя к непрерывному пределу $\textstyle \Delta t\to 0$, получаем ($\textstyle t_{0}=0$):

$g_{\mu \alpha }g_{\nu \alpha }=\sum _{i}G_{\mu \alpha }(t_{i})G_{\nu \alpha }(t_{i})\Delta t=\int \limits _{0}^{t}G_{\mu \alpha }(\tau )G_{\nu \alpha }(\tau )\,d\tau ,$

или:

$\mathbf {g} (t)\cdot \mathbf {g} ^{T}(t)=\int \limits _{0}^{t}e^{-\mathbf {A} \tau }\,\mathbf {B} \,\mathbf {B} ^{T}e^{-\mathbf {A} ^{T}\tau }\;d\tau .$

(6.26)

Напомню, что $\textstyle (\mathbf {A} \cdot \mathbf {B} )^{T}=\mathbf {B} ^{T}\cdot \mathbf {A} ^{T}$ (см. стр. \pageref{math_mat_tensor}). Решение для $\textstyle \mathbf {y}$ запишем в матричном виде, учитывая, что $\textstyle \mathbf {y} _{0}=\mathbf {x} _{0}$ при $\textstyle t=0$:

$\mathbf {y} =\mathbf {x} _{0}+\mathbf {g} (t)\cdot \mathbf {\epsilon }$

Поэтому, так как $\textstyle \mathbf {x} =e^{\mathbf {A} t}\,\mathbf {y}$, окончательное решение системы линейных стохастических уравнений имеет вид:

$\mathbf {x} (t)={\bar {\mathbf {x} }}(t)+\mathbf {S} (t)\cdot \mathbf {\epsilon } ,$

(6.27)

где $\textstyle \mathbf {S} =e^{\mathbf {A} t}\,\mathbf {g}$. Вектор $\textstyle \epsilon =\{\varepsilon _{1},...,\varepsilon _{n}\}$ представляет собой набор независимых случайных чисел с гауссовым распределением, имеющим нулевое среднее и единичную дисперсию, а $\textstyle {\bar {\mathbf {x} }}(t)$ — среднее значение (6.20), (6.24). В качестве упражнения ($\textstyle \lessdot$ H) предлагается найти матрицу $\textstyle e^{\mathbf {A} t}$ для двухмерного осциллятора и проверить решение (6.27).

$\textstyle \bullet$ Вычислим матрицу дисперсий:

$D_{\alpha \beta }=\left\langle (x-{\bar {x}})_{\alpha }(x-{\bar {x}})_{\beta }\right\rangle =S_{\alpha i}S_{\beta j}\left\langle \varepsilon _{i}\varepsilon _{j}\right\rangle =\left[\mathbf {S} \,\mathbf {S} ^{T}\right]_{\alpha \beta }=\left[e^{\mathbf {A} \,t}\,\mathbf {g} \,\mathbf {g} ^{T}\,e^{\mathbf {A} ^{T}\,t}\right]_{\alpha \beta }.$

Учитывая (6.26), имеем:

$\mathbf {D} (t)=\mathbf {S} \,\mathbf {S} ^{T}=\int \limits _{0}^{t}e^{\mathbf {A} (t-\tau )}\,\mathbf {B} \,\mathbf {B} ^{T}\,e^{\mathbf {A} ^{T}(t-\tau )}\,d\tau .$

(6.28)

Это соотношение можно ($\textstyle \lessdot$ H) сразу получить из уравнения для средних (6.17), из которых следует матричное уравнение:

${\dot {\mathbf {D} }}=\mathbf {A} \cdot \mathbf {D} +\mathbf {D} \cdot \mathbf {A} ^{T}+\,\mathbf {B} \cdot \mathbf {B} ^{T}.$

(6.29)

Если существует стационарный режим, то $\textstyle {\dot {\mathbf {D} }}=0$ и уравнение (6.29) позволяет легко найти $\textstyle \mathbf {D}$.

$\textstyle \bullet$ Распределение для $\textstyle x$ имеет гауссовый вид, поэтому, зная матрицу дисперсий, можно записать марковскую плотность вероятности:

$P(\mathbf {x} _{0},0\Rightarrow \mathbf {x} ,t)={\frac {1}{(2\pi )^{n/2}{\sqrt {\det \mathbf {D} (t)}}}}\exp \left[-{\frac {1}{2}}(x-{\bar {x}})_{\alpha }\,D_{\alpha \beta }^{-1}(t)\,(x-{\bar {x}})_{\beta }\right],$

где $\textstyle \mathbf {D} ^{-1}$ — обратная матрица дисперсий и $\textstyle {\bar {\mathbf {x} }}={\bar {\mathbf {x} }}(t)$ — средние значения динамических переменных. Они полностью определяют свойства процесса. В частности, характеристическая функция ($\textstyle \lessdot$ H):

$\left\langle e^{\imath \,\mathbf {p} \cdot \mathbf {x} }\right\rangle =e^{\imath \,\mathbf {p} \cdot {\bar {\mathbf {x} }}-{\frac {1}{2}}\,\mathbf {p} \cdot \mathbf {D} \cdot \mathbf {p} }$

позволяет легко находить моменты произвольных порядков.

$\textstyle \bullet$ При помощи (6.27), (6.28) несложно ($\textstyle \lessdot$ H) найти ковариационную матрицу:

$\mathrm {cov} \,{\alpha \beta }(t,t+\tau )=\left\langle x_{\alpha }(t)x_{\beta }(t+\tau )\right\rangle -\left\langle x_{\alpha }(t)\right\rangle \left\langle x_{\beta }(t+\tau )\right\rangle =\mathbf {D} (t)\,e^{\mathbf {A} ^{T}\,\tau }.$

(6.30)

Если в пределе $\textstyle t\to \infty$ у системы существует стационарный режим, то в этом случае матрица дисперсий $\textstyle \mathbf {D}$ становится постоянной, а ковариация зависит только от разности времён $\textstyle \tau$.

$\textstyle \bullet$ Таким образом, алгоритм решения линейной задачи следующий:

• $\textstyle \triangleright$ Находим собственные значения и вектора матрицы $\textstyle \mathbf {A}$.

• $\textstyle \triangleright$ Записываем решение для средних (6.24) и выражаем $\textstyle \mu _{k}$ через $\textstyle \mathbf {x} _{0}$.

• $\textstyle \triangleright$ При помощи соотношения $\textstyle \left[e^{\mathbf {A} t}\right]_{\alpha \beta }=\partial {\overline {x}}_{\alpha }/x_{0\beta }$ находим $\textstyle e^{\mathbf {A} t}$.

• $\textstyle \triangleright$ Вычисляем матрицу дисперсий $\textstyle D_{\alpha \beta }$.

---

Уравнение стохастического осциллятора <<	Оглавление	>> Многомерие помогает одномерию

---

Стохастический мир - простое введение в стохастические дифференциальные уравнения