Восстановление итерационной схемы
по плотности инвариантной меры

Ограничимся динамическими процессами в виде итерационных схем

$$x_{n+1}=f(x_n)\,,n=0,1,\ldots,$$ (1)

в которых $n$ играет роль времени, $x_n\in[0,1], f(x)\in[0,1], f(0)=f(1)=0.$ Мы уже знаем из § 1, что если $f(x)$ непрерывная в $[0,1]$ функция от$x$ , имеющая единственный максимум в $(0,1)$ , то, как правило, при увеличении $n$ поведение $x_n$ становится непредсказуемым, хаотическим: последовательные итерации (1) распределяются в $[0,1]$ случайно с некоторой плотностью вероятности $\rho (x)$ , не зависящей ни от $n$ , ни от начального значения $x_0$ (и потому называемой плотностью инвариантной меры). Данное свойство динамической системы (1) используется для машинной реализации датчика случайных чисел с плотностью распределения $\rho (x)$ в тех (весьма немногочисленных) случаях,когда $\rho (x)$ является плотностью инвариантной меры известного динамического закона $f(x)$ . Поэтому приобретает большое значение задача отыскания итерационной схемы (1), для которой заданная $\rho (x)$ служит плотностью инвариантной меры, называемая нами простейшей обратной задачей. Соответствующая прямая задача, когда по известной $f(x)$ определяется $\rho (x)$ , решается уравнением Фробениуса - Перрона (1.10)

$$\rho (x)=\int_0^1\,\rho (y)\,\delta (x-f(y))\,dy\,,$$ (2)

при выводе которого используется марковское свойство итерационной процедуры (1)с параметром $n$ как временем и независимость $\rho (x)$ от $n$ и $x_0$ .

Уравнение (2) решает и обратную задачу, если поменять в нем местами известную и неизвестную функции. Благодаря наличию $\delta$ -функции Дирака под интегралом уравнение (2) можно превратить в функциональное (для прямой задачи) или дифференциальное (для обратной) следующим образом. Воспользуемся одновершинностьюфункции $f(x)$ и представим ее в виде

$$f(x)=\cases {l(x)\quad\textrm{при}\quad 0\leq x\leq x_m\,,\cr p(1-x)\quad\textrm{при}\quad x_m\leq x\leq 1\,,\cr}$$ (3)

введя левую $l(x)$ и правую $p(1-x)$ ветви функции $f(x)$ .Здесь $l(x)$ , $p(x)$ -- монотонно возрастающие функции, $x_m$ -- точка максимума $f(x)$ . Очевидно,

$$l(x_m)=p(1-x_m),\qquad l(0)=p(0)=0\,,$$ (4)

поскольку $f(x)$ непрерывна и на концах интервала $[0,1]$ обращается в нуль.

Подставив (3) в (2) и воспользовавшись (4), получаем уравнение Фробениуса - Перрона в виде

$$\rho (x)=\int_0^{x_m}[\rho (y)\delta (x-l(y))+\rho (1-y)\delta (x-p(y))]\,dy\,,$$

откуда, ликвидировав интегрирование с помощью $\delta$ -функции по формуле (1.9), находим

$$\rho (x)={\frac{\rho (l^{-1}(x))}{l'(l^{-1}(x))}}+{\frac{\rho... ...))}{p'(p^{-1}(x))}}\quad\textrm{при}\quad 0\leq x\leq l(x_m)\,,$$ (5)

где $l^{-1}(x)$ ,$p^{-1}(x)$ -- (очевидно, однозначные) функции, обратные$l(x)$ , $p(x)$ соответственно, а штрихом отмечается дифференцирование по аргументу.

Равенство (5) при известной $f(x)$ , т.е. известных $l(x)$ , $p(x)$ и $x_m$ ,-- функциональное уравнение относительно $\rho (x)$ . Поэтому столь труднорешаема прямая задача (см. примеры в § 1). Для обратной же задачи, когда известна $\rho (x)$ , а неизвестны $l(x)$ , $p(x)$ и $x_m\in (0,1)$ , равенство (5) представляет собой дифференциальное уравнение первого порядка, решение которого не является проблемой. Если, вдобавок, $l(x)$ и $p(x)$ связаны известной аналитической зависимостью, то решение может быть и единственным.

В качестве дополнительных соображений, могущих обеспечить указанную аналитическую связь, простейшими, видимо, являются свойства симметрии (или асимметрии) функции $f(x)$ , задаваемые априори. Мы ограничимся следующим представлением:

$$p(x)=l(g(x))\,,$$ (6)

где $g(x)$ -- известная дифференцируемая монотонная функция, $g(0)=0$ .Уравнение (6) назовем уравнением асимметрии. Оно играет здесь ту же роль, какую играет уравнение состояния в механике сплошных сред, -- роль замыкающего соотношения, обеспечивающего совпадение числа независимых уравнений с числом неизвестных функций.

В силу однозначности $l(x)$ из (4) и (6) получаем

$$x_m=g(1-x_m)\,,$$ (7)

откуда определяется $x_m$ . Используя (6), приводим (5) к виду

$$l'(l^{-1}(x))\rho (x)=\rho (l^{-1}(x))+{\frac{\rho (1-g^{-1}... ...-1}(x)))}{g'(g^{-1}(l^{-1}(x)))}}\,,\quad 0\leq x\leq l(x_m)\,,$$

где $g^{-1}(x)$ -- функция, обратна$g(x)$ (она однозначна). Если в последнем равенстве ввести новую переменную $y=l^{-1}(x)$ , то получим обыкновенное дифференциальное уравнение первого порядка относительно $l(y)$ с разделяющимися переменными:

$$\rho (l(y))l'(y)=\rho (y)+\rho (1-g^{-1}(y))/g'(g^{-1}(y))\,,\qquad 0\leq y\leq x_m\,.$$ (8)

Его нужно решать при начальном условии $l(0)=0.$ Когда $\rho (x)$ и $g(x)$ известны, решение (8) в интервале $y\in [0,x_m]$ представляется без труда[9, с.35]. Для определения правой ветви функции$f(x)$ используется (7). Мы выпишем решение в некоторых интересных частных случаях $1^0$ -$3^0$ ипроиллюстрируем примерами.

Ясно из (8), что инвариантную плотность $\rho (x)$ достаточно задать с точностью до постоянного множителя (обычно его выбирают из условия нормировки (1.22)).

$1^0.$ Функция $f(x)$ симметрична относительно точки $x=1/2$ , т.е.$f(x)=f(1-x).$ В этом случае (6) принимает вид $p(x)=l(x)$ , так что$g(x)=g^{-1}(x)=x,x_m=1/2,$ и мы имеем из равенства (8), заменив в нем $y$ на$x$ (т.е. переобозначив аргумент),

$$\rho (l(x))l'(x)=\rho (x)+\rho (1-x)\,,\qquad 0\leq x\leq 1/2\,.$$

Если, кроме того, $\rho (x)=\rho (1-x)$ , то (8) упростится еще более:

$$\rho (l(x))l'(x)=2\rho (x)\,,\qquad 0\leq x\leq 1/2\,.$$ (9)

П р и м е р 1. Пусть $\rho (x)=x^\alpha$ при $0\leq x\leq 1/2$ , $\rho (x)=(1-x)^\alpha$ при$1/2\leq x\leq 1\,,\alpha >-1.$ Так как $\rho (x)=\rho (1-x)$ , естественно предположить,что и $f(x)=f(1-x)$ . Интегрируя (9), находим $l(x)=2^{1/(1+\alpha )}x$ при$0\leq x\leq 1/2$ и в силу уравнения асимметрии (6)

$$f(x)=2^{\frac{1}{\alpha +1}}\cases {x\,,&$\qquad 0\leq x\leq 1/2\,,$\cr 1-x\,,&$\qquad 1/2\leq x\leq 1\,.$\cr}$$

Так как$f(x)\leq 1$ , необходимо $\alpha \geq 0$ .

Этот пример интересен тем, что определяет плотность инвариантной меры хаотического аттрактора треугольного отображения (1.17)

$$x_{n+1}=2r\cdot\cases {x_n\,,&$\qquad 0\leq x_n\leq 1/2\,,$\cr 1-x_n\,,&$\qquad 1/2\leq x_n\leq 1\,,$\cr}$$ (10)

при любом $r\in [1/2,1]$ , причем$\alpha$ находится из соотношения $2r=2^{1/(1+\alpha )}$ . В примере 2§ 1, решая прямую задачу, мы смогли найти $\rho (x)$ только при$r=1$ .

Треугольное отображение (10) демонстрирует жесткий переход к хаосу. Именно, при $0<r<1/2$ аттрактором является только точка $x=0$ ($\vert f'\vert=2r<1$ ), апри $1/2<r<1$ мы имеем уже хаотический аттрактор с плотностью инвариантной меры. Показатель Ляпунова для этого отображения вычисляется легко по формуле (1.6) и равен$\lambda =\ln 2r$ . Он отрицателен при $0<r<1/2$ и положителен при$1/2<r<1$ . Таким образом, первоначально близкие точки в итерационном процессе могут "экспоненциально разбежаться" с показателем экспоненты $n\ln 2r$ , т.е. первоначально спокойный режим $x=0$ при переходе через точку $r=1/2$ резко становится хаотическим, почему такого рода возникновение хаоса и называется жестким.

П р и м е р 2. Пусть $\rho (x)=x^{\alpha }(1-x)^{\alpha }$ , $\alpha >-1$ . Тогда$\rho (x)=\rho (1-x)$ и, снова предполагая, что $f(x)=f(1-x)$ , имеем из (9)

$$B_{l(x)}(\alpha +1,\alpha +1)=2B_x(\alpha +1,\alpha +1)\,,$$ (11)

где $B_x(\alpha ,\beta )$ -- неполная бета-функция [10]. При некоторых $\alpha$ уравнение (11) разрешимо в явных функциях относительно $l(x)$ . Так, если $\alpha =-1/2$ ,то из (11) получаем $f(x)=4x(1-x)$ , т.е. логистическое отображение (1.23).Если $\alpha =0,$ то вследствие (10) $l(x)=2x$ при $0\leq x\leq 1/2$ , так что всилу $f(x)=f(1-x)$ находим

$$f(x)=\cases {2x\,,&$\qquad 0\leq x\leq 1/2\,,$\cr 2(1-x)\,,&$\qquad 1/2\leq x\leq 1\,,$\cr}$$

т.е. треугольное отображение (10) при $r=1$ .

$2^0.$ Линейная асимметрия правой и левой ветвей $f(x)$ . Возьмем$g(x)$ в виде линейной функции: $g(x)=kx,k>0.$ Тогда в (7) $x_m$ также войдет линейно, и мы имеем $x_m=k/(k+1)$ ,

$$f(x)=\cases {l(x)\,,&$\qquad 0\leq x\leq x_m\,,$\cr l(k(1-x))\,,&$\qquad x_m\leq x\leq 1\,,$\cr}$$ (12)

т. е. уравнение асимметрии (6) примет вид $p(x)=l(kx).$ Так как$g_{-1}=x/k$ , $g'(x)=k,$ то уравнение (8) запишется в следующей форме($y\equiv x$ ):

$$\rho (l(x))l'(x)=\rho (x)+\rho (1-x/k)/k\,, \qquad 0\leq x\leq x_m\,.$$ (13)

При $k=1$ получаем отсюда предыдущий случай $1^0$ .

Когда $\rho (1-x/k)=\rho (x)$ , уравнение (13) упрощается:

$$x_m\rho (l(x))l'(x)=\rho (x)\,.$$ (14)

Пример 3. Пусть $\rho (x)=x^{\alpha }$ при $0\leq x\leq x_{\star}$ и$\rho (x)=k^{\alpha }(1-x)^{\alpha }$ при $x_{\star}\leq x\leq 1$ . Требуя непрерывность $\rho (x)$ , находим $x_{\star}=k/(1+k)$ , так что естественно предположить справедливость (12), причем $x_m=x_{\star}$ . Уравнение (14) примет вид $l^{\alpha }(x)l'(x)=(1+1/k)x^{\alpha }$ . Интегрируя его с начальным условием$l(0)=0$ ,получаем $l(x)=x_{\star}^{-1/(1+\alpha )}x,$ и в силу (12) функция $f(x)$ равна$x/x_{\star}^{1/(1+\alpha )}$ при $0\leq x\leq x_{\star}$ и равна$k(1-x)/x_{\star}^{1/(1+\alpha )}$ при $x_{\star}\leq x\leq 1$ .

Как и в примере 1, в этом примере с точностью до нормирующего множител определяется инвариантная плотность хаотического аттрактора динамическойсистемы

$$x_{n+1}=r\cases{ x_n\,,&$\qquad 0\leq x_n\leq x_{\star}\,,$\cr k(1-x_n)\,,&$\qquad x_{\star}\leq x_n\leq 1\,,$\cr}$$

при каждом $r\in [1,1+1/k],$ причем $\alpha$ определяется из уравнени$r^{1+\alpha }=1+1/k.$ Чтобы $f(x)\leq 1,$ необходимо $\alpha \geq 0\,.$

$3^0.$ Нелинейная асимметрия. Если в (8)

$$\rho (x)=\rho (1-g^{-1}(x))\,, \qquad 0\leq x\leq x_m\,,$$ (15)

где $g^{-1}(x)$ известная неотрицательная монотонная функция, причем$g^{-1}(0)=0$ ,то $x_m$ удовлетворяет соотношению $x_m=1-g^{-1}(x_m)$ . В уравнении асимметрии(6) в качестве $g(x)$ возьмем функцию, обратную $g^{-1}(x)$ . Равенство (8) поэтому запишется в виде

$$\rho (l(x))l'(x)=[1+1/g'(g^{-1}(x))]\rho (x)\,, \quad 0\leq x\leq x_m\,,$$ (16)

включающем в себя (9) и (14) частными случаями.

П р и м е р 4. Пусть $\rho (x)=x^{\alpha }$ при $0\leq x\leq x_m$ ,$\rho (x)=(1-x)^{\beta }$ при $x_m\leq x\leq 1$ , $\alpha \geq -1$ , $\beta \geq -1$ , $\alpha \beta >0$ . Из условия непрерывности$\rho (x)$ в точке $x_m$ получаем уравнение $x^{\alpha }_m+(1-x_m)^{\beta }=1$ для определения $x_m\in (0,1)$ . Возьмем $g^{-1}=y^{\alpha /\beta }$ , так что$g(y)=y^{\beta /\alpha }$ . Тогда справедливо (15), и из (16) следует

$$l^{\alpha }(x)l'(x)=(1+\alpha x^{\alpha /\beta -1}/\beta )x^{\alpha }\,, \qquad 0\leq x\leq x_m\,,$$

откуда определяется левая ветвь функции $f(x)$ :

$$l(x)=x\left(1+{\frac{1+\alpha }{1+\beta }}x^{\alpha /\beta -1}\right)^{1/(1+\alpha )}.$$

Правая ветвь восстанавливается по уравнению (6).

Рассмотренные примеры наталкивают на мысль, что угадать уравнение асимметрии можно по инвариантной плотности. Так, в примерах 1 и 2, пользуясь симметрией $\rho (x)$ относительно точки $x=1/2$ , мы предположили такую же симметрию для $f(x)$ , а в примере 3, опираясь на $\rho (x)=\rho (k(1-x))$ , $0\leq x\leq x_{\star}$ , записали уравнение асимметрии (6) в виде $p(x)=l(kx)$ . Тем не менее, ответ навопрос о том, единственно ли уравнение асимметрии для данной $\rho (x)$ , отрицателен. Это уже можно усмотреть в примерах 1 и 3 при $\alpha =0$ , когда для $\rho (x)=1$ имеем два разных уравнения асимметрии и, соответственно, два разных динамических закона$f(x)$ . Этот же ответ дает следующий нетривиальный пример.

П р и м е р 5. Возьмем $\rho (x)=x$ при $0\leq x\leq 1/2$ , $\rho (x)=1-x$ при$1/2\leq x\leq 1$ . Положим в (6) $g(x)=(x/a)^{1/\beta }$ , $x_m=1/2$ , $\beta >0$ ,$a>0$ . Тогда $g^{-1}(x)=ax^{\beta }$ , а из (7) находим $a=2^{\beta -1}$ . Интегрируя (8) и используя (6), получаем

$$f(x)=\cases{ [x^2+(a/\alpha )^2x^{2\beta }]^{1/2}\,,\ 0\leq x... ...^{-2/\beta }(1-x)^{2/\beta }]^{1/2}\,,\ 1/2\leq x\leq 1 \,.\cr}$$ (17)

В то же время, взяв $g(x)=x$ , в примере 1 для $\alpha =1$ имеем $f(x)=x\sqrt 2$ при $0\leq x\leq 1/2$ , $f(x)=(1-x)\sqrt 2$ при $1/2\leq x\leq 1$ , т. е. имеем функцию совершенно другого типа, чем (17).

Таким образом, единственность обратной задачи определяется уравнениемасимметрии (6): разные предположения о зависимости правой ветви от левойприводят к разным решениям обратной задачи, причем несимметричный динамическийзакон $f(x)$ может иметь симметричную инвариантную плотность.

З а м е ч а н и е. Покажем, как с помощью решения обратной задачи можно найтипримеры хаотизации известных вычислительных схем.

О хаотизации вычислительной схемы Ньютона при отыскании корней функции мы рассказали в § 1.

Более просто строится пример хаотизации в методе секущих решения уравнения $F(x)=0$ . Если итерации начинаются с точки $x_0$ , то итерационная схема имеетвид [11]

$$x_{n+1}=\phi (x_n,x_0),\qquad n=0,1,\ldots,$$

где

$$\phi(x_n,x_0)=[x_0F(x)-xF(x_0)]/[F(x)-F(x_0)]\,.$$ (18)

Приравняв (18) логистической функции $4rx(1-x)$ , получаем

$$F(x)=x(1-4r+4rx)/(x_0-4rx+4rx^2)\,,$$

так что корни уравнения $F(x)=0$ в $[0,1]$ суть $x=1-1/4r$ при $r\geq 1/4$ и $x=0$ ,а метод секущих с начальной точкой $x_0$ при $4r\geq 3$ снова приводит кпериодическим аттракторам или к хаотическому (посколькуэти аттракторынаблюдаются для логистического отображения).

Наконец, в основе решения задачи Коши

$$x'(t)=F(x(t)),\quad x(0)=x_0$$

методом Эйлера с постоянным шагом $h$ лежит итерационная процедура [11] $x_{n+1}=x_n+hF(x_n)$ , и нетрудно проверить, что при $F(x)=-x^2$ , $x(0)=0$ этапроцедура сводится к итерированию функции $y(1-y)/h$ , где $y=1-hx$ , т. е. китерированию логистического отображения. Как мы видели, уже при$h=1/4$ (даже немного ранее) появляется хаотический аттрактор, так что к решению $x(t)=x_0/(1+tx_0)$ задачи эта схема не сходится.