2014 dxdy logo

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

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


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


Посмотреть правила форума



Начать новую тему Ответить на тему На страницу 1, 2, 3  След.
 
 численные методы для траекторий хаотических систем
Сообщение27.05.2017, 18:22 


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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение27.05.2017, 23:45 
Заслуженный участник


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

(Оффтоп)

код: [ скачать ] [ спрятать ]
Используется синтаксис Matlab M
function lorenz(action)
%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 (hence the word "strange").
%
%   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;
 

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение28.05.2017, 00:44 


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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение28.05.2017, 10:18 
Заслуженный участник


25/02/11
1786
Igor_Dmitriev в сообщении #1219164 писал(а):
В то же время, по-моему, есть методы, позволяющие этого избежать.

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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение28.05.2017, 10:28 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение28.05.2017, 14:28 
Заслуженный участник


05/08/14
1564
Igor_Dmitriev в сообщении #1219283 писал(а):
Только вот при изменении начальных условий на сотые доли процента траектория сильно меняется.

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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение28.05.2017, 15:30 


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

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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение28.05.2017, 17:53 
Заслуженный участник


05/08/14
1564
Igor_Dmitriev в сообщении #1219407 писал(а):
А как тогда решили вопрос с состоянием солнечной системы на миллиарды лет в прошлом и будущем?
Взяли суперкомпьютер, выбрали алгоритм, контролирующий точность на каждом шаге, запустили программу, через месяц получили решение.

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение01.06.2017, 03:01 


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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение01.06.2017, 04:07 
Аватара пользователя


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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение01.06.2017, 04:39 


15/01/12
196
$x$ я получил в матлаб из $ode45$, а потом сделал
Код:
plot(filter(1, [1 -1], x-mean(x)))

Изображение

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение01.06.2017, 07:29 
Заслуженный участник


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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение01.06.2017, 15:48 


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

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение01.06.2017, 16:06 
Заслуженный участник


09/05/12
25179
Igor_Dmitriev в сообщении #1221006 писал(а):
На той картинке, которую я привел, периодов примерно $10^4-10^5$.
А пикселей по горизонтали у монитора $\sim 10^3$. Естественно, что "гладкость" картинки при этом никакого отношения к реальности не имеет.

 Профиль  
                  
 
 Re: численные методы для траекторий хаотических систем
Сообщение01.06.2017, 16:30 


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 36 ]  На страницу 1, 2, 3  След.

Модераторы: Модераторы Математики, Супермодераторы



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

Сейчас этот форум просматривают: Vladimir Pliassov


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group