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  След.

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



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

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


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

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