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

Материал из Synset

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

$\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}$ можно разбить на сумму двух матриц:

$\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^{(k)}_\alpha \, u^{*(k)}_\beta \,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^t_{0} G_{\mu\alpha}(\tau)G_{\nu\alpha}(\tau)\,d\tau,$

или:

$\mathbf{g}(t)\cdot \mathbf{g}^T(t) = \int\limits^t_{0} 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^t_0 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^{-1}_{\alpha\beta}(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}$ .

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

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

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

Материал из Synset

Просмотры

Личные инструменты

Навигация

почитай

Поиск

Инструменты