численные методы для траекторий хаотических систем

Igor_Dmitriev · 15/01/12 215

Как известно, малые отклонения, произошедшие в достаточно удалённом моменте времени, приводят в случае хаотических систем к существенным отклонениям.
Кроме того, погрешности возникают хотя бы по причине конечной точности на компьютере.
В то же время, по-моему, есть методы, позволяющие этого избежать.
Что это за методы?

dsge · 05/08/14 1564

Посмотрите какой-нибудь пример, например, в Матлабе

(Оффтоп)

код: [ скачать ] [ спрятать ]

Используется синтаксис Matlab M

function lorenz&#40;action&#41;

%LORENZ Plot the orbit around the Lorenz chaotic attractor.

%   This demo animates the integration of  the

%   three coupled nonlinear differential equations

%   that define the "Lorenz Attractor", a chaotic

%   system first described by Edward Lorenz of

%   the Massachusetts Institute of Technology.

%

%   As the integration proceeds you will see a

%   point moving around in a curious orbit in

%   3-D space known as a strange attractor. The

%   orbit is bounded, but not periodic and not

%   convergent &#40;hence the word "strange"&#41;.

%

%   Use the "Start" and "Stop" buttons to control

%   the animation.

%   Adapted for Demo by Ned Gulley, 6-21-93; jae Roh, 10-15-96

%   Copyright 1984-2005 The MathWorks, Inc.

%   $Revision: 5.13.4.3 $  $Date: 2005/12/15 20:52:53 $

% The values of the global parameters are

global SIGMA RHO BETA

SIGMA = 10.;

RHO = 28.;

BETA = 8./3.;

% Possible actions:

% initialize

% close

% Information regarding the play status will be held in

% the axis user data according to the following table:

play= 1;

if nargin<1,

    action=&#39;initialize&#39;;

end

switch action

    case &#39;initialize&#39;

        oldFigNumber=watchon;

        figNumber=figure&#40; ...

            &#39;Name&#39;,&#39;Lorenz Attractor&#39;, ...

            &#39;NumberTitle&#39;,&#39;off&#39;, ...

            &#39;Visible&#39;,&#39;off&#39;&#41;;

        colordef&#40;figNumber,&#39;black&#39;&#41;

        axes&#40; ...

            &#39;Units&#39;,&#39;normalized&#39;, ...

            &#39;Position&#39;,[0.05 0.10 0.75 0.95], ...

            &#39;Visible&#39;,&#39;off&#39;&#41;;

        text&#40;0,0,&#39;Press the "Start" button to see the Lorenz demo&#39;, ...

            &#39;HorizontalAlignment&#39;,&#39;center&#39;&#41;;

        axis&#40;[-1 1 -1 1]&#41;;

        %===================================

        % Information for all buttons

        xPos=0.85;

        btnLen=0.10;

        btnWid=0.10;

        % Spacing between the button and the next command&#39;s label

        spacing=0.05;

        %====================================

        % The CONSOLE frame

        frmBorder=0.02;

        yPos=0.05-frmBorder;

        frmPos=[xPos-frmBorder yPos btnLen+2*frmBorder 0.9+2*frmBorder];

        uicontrol&#40; ...

            &#39;Style&#39;,&#39;frame&#39;, ...

            &#39;Units&#39;,&#39;normalized&#39;, ...

            &#39;Position&#39;,frmPos, ...

            &#39;BackgroundColor&#39;,[0.50 0.50 0.50]&#41;;

        %====================================

        % The START button

        btnNumber=1;

        yPos=0.90-&#40;btnNumber-1&#41;*&#40;btnWid+spacing&#41;;

        labelStr=&#39;Start&#39;;

        callbackStr=&#39;lorenz&#40;&#39;&#39;start&#39;&#39;&#41;;&#39;;

        % Generic button information

        btnPos=[xPos yPos-spacing btnLen btnWid];

        startHndl=uicontrol&#40; ...

            &#39;Style&#39;,&#39;pushbutton&#39;, ...

            &#39;Units&#39;,&#39;normalized&#39;, ...

            &#39;Position&#39;,btnPos, ...

            &#39;String&#39;,labelStr, ...

            &#39;Interruptible&#39;,&#39;on&#39;, ...

            &#39;Callback&#39;,callbackStr&#41;;

        %====================================

        % The STOP button

        btnNumber=2;

        yPos=0.90-&#40;btnNumber-1&#41;*&#40;btnWid+spacing&#41;;

        labelStr=&#39;Stop&#39;;

        % Setting userdata to -1 &#40;=stop&#41; will stop the demo.

        callbackStr=&#39;set&#40;gca,&#39;&#39;Userdata&#39;&#39;,-1&#41;&#39;;

        % Generic  button information

        btnPos=[xPos yPos-spacing btnLen btnWid];

        stopHndl=uicontrol&#40; ...

            &#39;Style&#39;,&#39;pushbutton&#39;, ...

            &#39;Units&#39;,&#39;normalized&#39;, ...

            &#39;Position&#39;,btnPos, ...

            &#39;Enable&#39;,&#39;off&#39;, ...

            &#39;String&#39;,labelStr, ...

            &#39;Callback&#39;,callbackStr&#41;;

        %====================================

        % The INFO button

        labelStr=&#39;Info&#39;;

        callbackStr=&#39;lorenz&#40;&#39;&#39;info&#39;&#39;&#41;&#39;;

        infoHndl=uicontrol&#40; ...

            &#39;Style&#39;,&#39;push&#39;, ...

            &#39;Units&#39;,&#39;normalized&#39;, ...

            &#39;position&#39;,[xPos 0.20 btnLen 0.10], ...

            &#39;string&#39;,labelStr, ...

            &#39;call&#39;,callbackStr&#41;;

        %====================================

        % The CLOSE button

        labelStr=&#39;Close&#39;;

        callbackStr= &#39;close&#40;gcf&#41;&#39;;

        closeHndl=uicontrol&#40; ...

            &#39;Style&#39;,&#39;push&#39;, ...

            &#39;Units&#39;,&#39;normalized&#39;, ...

            &#39;position&#39;,[xPos 0.05 btnLen 0.10], ...

            &#39;string&#39;,labelStr, ...

            &#39;call&#39;,callbackStr&#41;;

        % Uncover the figure

        hndlList=[startHndl stopHndl infoHndl closeHndl];

        set&#40;figNumber,&#39;Visible&#39;,&#39;on&#39;, ...

            &#39;UserData&#39;,hndlList&#41;;

        set&#40;figNumber, &#39;CloseRequestFcn&#39;, &#39;clear global SIGMA RHO BETA;closereq&#39;&#41;;

        watchoff&#40;oldFigNumber&#41;;

        figure&#40;figNumber&#41;;

    case &#39;start&#39;

        axHndl=gca;

        figNumber=gcf;

        hndlList=get&#40;figNumber,&#39;UserData&#39;&#41;;

        startHndl=hndlList&#40;1&#41;;

        stopHndl=hndlList&#40;2&#41;;

        infoHndl=hndlList&#40;3&#41;;

        closeHndl=hndlList&#40;4&#41;;

        set&#40;[startHndl closeHndl infoHndl],&#39;Enable&#39;,&#39;off&#39;&#41;;

        set&#40;stopHndl,&#39;Enable&#39;,&#39;on&#39;&#41;;

        % ====== Start of Demo

        set&#40;figNumber,&#39;Backingstore&#39;,&#39;off&#39;&#41;;

        % The graphics axis limits are set to values known

        % to contain the solution.

        set&#40;axHndl, ...

            &#39;XLim&#39;,[0 40],&#39;YLim&#39;,[-35 10],&#39;ZLim&#39;,[-10 40], ...

            &#39;Userdata&#39;,play, ...

            &#39;XTick&#39;,[],&#39;YTick&#39;,[],&#39;ZTick&#39;,[], ...

            &#39;Drawmode&#39;,&#39;fast&#39;, ...

            &#39;Visible&#39;,&#39;on&#39;, ...

            &#39;NextPlot&#39;,&#39;add&#39;, ...

            &#39;Userdata&#39;,play, ...

            &#39;View&#39;,[-37.5,30]&#41;;

        xlabel&#40;&#39;X&#39;&#41;;

        ylabel&#40;&#39;Y&#39;&#41;;

        zlabel&#40;&#39;Z&#39;&#41;;

        % The orbit ranges chaotically back and forth around two different points,

        % or attractors.  It is bounded, but not periodic and not convergent.

        % The numerical integration, and the display of the evolving solution,

        % are handled by the function ODE23P.

        FunFcn=&#39;lorenzeq&#39;;

        % The initial conditions below will produce good results

        %y0 = [20 5 -5];

        % Random initial conditions

        y0&#40;1&#41;=rand*30+5;

        y0&#40;2&#41;=rand*35-30;

        y0&#40;3&#41;=rand*40-5;

        t0=0;

        tfinal=100;

        pow = 1/3;

        tol = 0.001;

        t = t0;

        hmax = &#40;tfinal - t&#41;/5;

        hmin = &#40;tfinal - t&#41;/200000;

        h = &#40;tfinal - t&#41;/100;

        y = y0&#40;:&#41;;

        % Save L steps and plot like a comet tail.

        L = 50;

        Y = y*ones&#40;1,L&#41;;

        cla;

        head = line&#40; ...

            &#39;color&#39;,&#39;r&#39;, ...

            &#39;Marker&#39;,&#39;.&#39;, ...

            &#39;markersize&#39;,25, ...

            &#39;erase&#39;,&#39;xor&#39;, ...

            &#39;xdata&#39;,y&#40;1&#41;,&#39;ydata&#39;,y&#40;2&#41;,&#39;zdata&#39;,y&#40;3&#41;&#41;;

        body = line&#40; ...

            &#39;color&#39;,&#39;y&#39;, ...

            &#39;LineStyle&#39;,&#39;-&#39;, ...

            &#39;erase&#39;,&#39;none&#39;, ...

            &#39;xdata&#39;,[],&#39;ydata&#39;,[],&#39;zdata&#39;,[]&#41;;

        tail=line&#40; ...

            &#39;color&#39;,&#39;b&#39;, ...

            &#39;LineStyle&#39;,&#39;-&#39;, ...

            &#39;erase&#39;,&#39;none&#39;, ...

            &#39;xdata&#39;,[],&#39;ydata&#39;,[],&#39;zdata&#39;,[]&#41;;

        % The main loop

        while &#40;get&#40;axHndl,&#39;Userdata&#39;&#41;==play&#41; && &#40;h >= hmin&#41;

            if t + h > tfinal, h = tfinal - t; end

            % Compute the slopes

            s1 = feval&#40;FunFcn, t, y&#41;;

            s2 = feval&#40;FunFcn, t+h, y+h*s1&#41;;

            s3 = feval&#40;FunFcn, t+h/2, y+h*&#40;s1+s2&#41;/4&#41;;

            % Estimate the error and the acceptable error

            delta = norm&#40;h*&#40;s1 - 2*s3 + s2&#41;/3,&#39;inf&#39;&#41;;

            tau = tol*max&#40;norm&#40;y,&#39;inf&#39;&#41;,1.0&#41;;

            % Update the solution only if the error is acceptable

            if delta <= tau

                t = t + h;

                y = y + h*&#40;s1 + 4*s3 + s2&#41;/6;

                % Update the plot

                Y = [y Y&#40;:,1:L-1&#41;];

                set&#40;head,&#39;xdata&#39;,Y&#40;1,1&#41;,&#39;ydata&#39;,Y&#40;2,1&#41;,&#39;zdata&#39;,Y&#40;3,1&#41;&#41;

                set&#40;body,&#39;xdata&#39;,Y&#40;1,1:2&#41;,&#39;ydata&#39;,Y&#40;2,1:2&#41;,&#39;zdata&#39;,Y&#40;3,1:2&#41;&#41;

                set&#40;tail,&#39;xdata&#39;,Y&#40;1,L-1:L&#41;,&#39;ydata&#39;,Y&#40;2,L-1:L&#41;,&#39;zdata&#39;,Y&#40;3,L-1:L&#41;&#41;

                drawnow;

            end

            % Update the step size

            if delta ~= 0.0

                h = min&#40;hmax, 0.9*h*&#40;tau/delta&#41;^pow&#41;;

            end

            % Bail out if the figure was closed.

            if ~ishandle&#40;axHndl&#41;

                return

            end

        end    % Main loop ...

        % ====== End of Demo

        set&#40;[startHndl closeHndl infoHndl],&#39;Enable&#39;,&#39;on&#39;&#41;;

        set&#40;stopHndl,&#39;Enable&#39;,&#39;off&#39;&#41;;

    case &#39;info&#39;

        helpwin&#40;mfilename&#41;;

end    % if strcmp&#40;action, ...

%===============================================================================

function ydot = lorenzeq&#40;t,y&#41;

%LORENZEQ Equation of the Lorenz chaotic attractor.

%   ydot = lorenzeq&#40;t,y&#41;.

%   The differential equation is written in almost linear form.

global SIGMA RHO BETA

A = [ -BETA    0     y&#40;2&#41;

    0  -SIGMA   SIGMA

    -y&#40;2&#41;   RHO    -1  ];

ydot = A*y;

Igor_Dmitriev · 15/01/12 215

Примеры смотрел, в той же википедии система Лоренца в матлабе решается с помощью метода Рунге-Кутта 4-5 порядка ode45. Только вот при изменении начальных условий на сотые доли процента траектория сильно меняется.

Vince Diesel · 25/02/11 1800

Igor_Dmitriev в сообщении #1219164 писал(а):

В то же время, по-моему, есть методы, позволяющие этого избежать.

Для фиксированного промежутка времени, если производить вычисления с достаточно большим количеством знаков после запятой, то такой проблемы (скажем, в системе Лоренца) не будет.

Евгений Машеров · 11/03/08 10041 Москва

Я бы начал с того, что хаотические системы это такие системы, в которых малые начальные возмущения будут приводить к большим отклонениям всегда, в силу их природы, а не несовершенства методов расчёта.
То есть всего, чего можно добиться - это получить зависимость результата только от возмущения начальных значений, а не ошибок вычисления. При этом ошибки вычисления на данный момент всегда можно рассматривать, как возмущения начальных значений при расчёте, начинающемся в данный момент. То есть нужно вести вычисления без ошибок. Если модель ограничена четырьмя действиями арифметики, то можно вести расчёты в рациональных числах, но числитель и знаменатель будут расти, и потребуется арифметика произвольной точности, которая рано или поздно исчерпает доступную память.
Поскольку начальные условия всегда задаются с конечной точностью (по крайней мере в прикладных задачах, где они получаются измерениями), смысл такого усложнения не вполне ясен.

dsge · 05/08/14 1564

Igor_Dmitriev в сообщении #1219283 писал(а):

Только вот при изменении начальных условий на сотые доли процента траектория сильно меняется.

Траектория начинает сильно (заметно) меняться через времена порядка $-\frac{\ln(0.0001)}{\lambda }$ , где $\lambda$ - старший показатель Ляпунова, т.е. это может быть довольно длительный период. Выбором шага разностной схемы (возможно, меняющемся, т.е. адаптированном к тому, где находится решение) и точности вычисления всегда можно добиться любой наперед заданной точности решения на любом большом, но фиксированном, интервале времени.

Igor_Dmitriev · 15/01/12 215

Спасибо за ответы.
Если я правильно понял -- то ничего не поделаешь: либо работать на конечном, заранее заданном промежутке, либо, при наличии лишь арифметических операций, иметь дело с рациональными выражениями, которые тоже ограничены ввиду быстрого роста знаменателя.

А как тогда решили вопрос с состоянием солнечной системы на миллиарды лет в прошлом и будущем?

dsge · 05/08/14 1564

Igor_Dmitriev в сообщении #1219407 писал(а):

А как тогда решили вопрос с состоянием солнечной системы на миллиарды лет в прошлом и будущем?

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

Igor_Dmitriev · 15/01/12 215

Чтобы не плодить темы, спрошу здесь же следующее:
есть уравнение, похожее на осцилятор Дуффинга
$\ddot{x} = x-x^3$
Если рассмотреть величину $\int\limits_{0}^{y}xdx$ за достаточно большой промежуток времени, то можно увидеть, что траектория $y$ ломанная, однако размерность Хаусдорфа близка к 0.
Почему интеграл траектории решения уравнения $\ddot{x} = x-x^3$ ломанная?

Aritaborian · 11/06/12 10390 стихия.вздох.мюсли

А где вы там ломаную увидели? Без картинки не совсем ясно, что у вас к чему. Может, это у вас точность вычислений хромает.

Igor_Dmitriev · 15/01/12 215

$x$ я получил в матлаб из $ode45$ , а потом сделал

Код:

plot(filter(1, [1 -1], x-mean(x)))

dsge · 05/08/14 1564

Такого не может быть, график должет быть гладкой периодической функцией или типа периодической с трендом. Сравните ваше численное решение с аналитическим - эллиптической функцией Якоби.
Интеграл у вас по $dt$ или $dx$ ? В последнем случае получим $y^2/2$ , что тоже должно быть гладко-периодично.

Igor_Dmitriev · 15/01/12 215

Интеграл по $dt$ , опечатка.
На малых интервалах, где бы они не были взяты, график действительно выглядит в виде гладкой периодической функции, при некоторых параметрах график также преобразуется в периодическую с трендом. На той картинке, которую я привел, периодов примерно $10^4-10^5$ .

Pphantom · 09/05/12 25179

Igor_Dmitriev в сообщении #1221006 писал(а):

На той картинке, которую я привел, периодов примерно $10^4-10^5$ .

А пикселей по горизонтали у монитора $\sim 10^3$ . Естественно, что "гладкость" картинки при этом никакого отношения к реальности не имеет.

Igor_Dmitriev · 15/01/12 215

Я знаю, что монитора недостаточно для полного отображения.
Вопрос в том, почему картинка похожа на броуновское движение.

Научный форум dxdy

Правила форума

численные методы для траекторий хаотических систем

Кто сейчас на конференции