Система Лоренца: определение и числовое решение

Система Лоренца - система обыкновенных дифференциальных уравнений (см. система Лоренца). Для реальных констант $$\sigma \mathrm{,\rho }\mathrm{,\beta }$$, система

Ценности Лоренца параметров для чувствительной системы - $$\sigma =\mathrm{}10,\beta =\mathrm{}\mathrm{8/3},\rho =\mathrm{}28$$. Запуск системы с [x(0),y(0),z(0)] = [10,20,10] и просмотр эволюции системы от времени 0 до 100.

sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];
xt0 = [10,20,10];
[tspan,a] = ode45(f,[0 100],xt0);     % Runge-Kutta 4th/5th order ODE solver
figure
plot3(a(:,1),a(:,2),a(:,3))
view([-10.0 -2.0])

Эволюция довольно сложная. Но за небольшой промежуток времени это выглядит как равномерное круговое движение. Постройте график решения по интервалу времени [0,1/10].

[tspan,a] = ode45(f,[0 1/10],xt0);     % Runge-Kutta 4th/5th order ODE solver
figure
plot3(a(:,1),a(:,2),a(:,3))
view([-30 -70])

Подбор кругового пути к решению ОДУ

Уравнения кругового пути имеют несколько параметров:

По этим параметрам определите положение круговой траектории на время xdata.

type fitlorenzfn

function f = fitlorenzfn(x,xdata)

theta = x(1:2);
R = x(3);
V = x(4);
t0 = x(5);
delta = x(6:8);
f = zeros(length(xdata),3);
f(:,3) = R*sin(theta(1))*sin(V*(xdata - t0)) + delta(3);
f(:,1) = R*cos(V*(xdata - t0))*cos(theta(2)) ...
    - R*sin(V*(xdata - t0))*cos(theta(1))*sin(theta(2)) + delta(1);
f(:,2) = R*sin(V*(xdata - t0))*cos(theta(1))*cos(theta(2)) ...
    - R*cos(V*(xdata - t0))*sin(theta(2)) + delta(2);

Чтобы найти наиболее подходящий круговой путь к системе Лоренца в моменты времени, указанные в решении ОДУ, используйте lsqcurvefit. Чтобы сохранить параметры в разумных пределах, установите ограничения на различные параметры.

lb = [-pi/2,-pi,5,-15,-pi,-40,-40,-40];
ub = [pi/2,pi,60,15,pi,40,40,40];
theta0 = [0;0];
R0 = 20;
V0 = 1;
t0 = 0;
delta0 = zeros(3,1);
x0 = [theta0;R0;V0;t0;delta0];
[xbest,resnorm,residual] = lsqcurvefit(@fitlorenzfn,x0,tspan,a,lb,ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

Постройте график наиболее подходящих круговых точек в момент времени из решения ОДУ вместе с решением системы Лоренца.

soln = a + residual;
hold on
plot3(soln(:,1),soln(:,2),soln(:,3),'r')
legend('ODE Solution','Circular Arc')
hold off

figure
plot3(a(:,1),a(:,2),a(:,3),'b.','MarkerSize',10)
hold on
plot3(soln(:,1),soln(:,2),soln(:,3),'rx','MarkerSize',10)
legend('ODE Solution','Circular Arc')
hold off

Подогнать ОДУ к дуге окружности

Теперь внесите изменения в параметры $$\lambda ,\beta$$и start, чтобы наилучшим образом соответствовать дуге окружности. Для еще лучшей подгонки позвольте изменить и начальную точку [10,20,10].

Для этого необходимо записать файл функции. paramfun который принимает параметры посадки ОДУ и вычисляет траекторию по времени t.

type paramfun

function pos = paramfun(x,tspan)

sigma = x(1);
beta = x(2);
rho = x(3);
xt0 = x(4:6);
f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];
[~,pos] = ode45(f,tspan,xt0);

Чтобы найти наилучшие параметры, используйте lsqcurvefit минимизировать различия между новой рассчитанной траекторией ОДУ и дугой окружности soln.

xt0 = zeros(1,6);
xt0(1) = sigma;
xt0(2) = beta;
xt0(3) = rho;
xt0(4:6) = soln(1,:);
[pbest,presnorm,presidual,exitflag,output] = lsqcurvefit(@paramfun,xt0,tspan,soln);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

Определите, насколько эта оптимизация изменила параметры.

fprintf('New parameters: %f, %f, %f',pbest(1:3))

New parameters: 9.132446, 2.854998, 27.937986

fprintf('Original parameters: %f, %f, %f',[sigma,beta,rho])

Original parameters: 10.000000, 2.666667, 28.000000

Параметры sigma и beta изменился примерно на 10%.

Постройте график измененного решения.

figure
hold on
odesl = presidual + soln;
plot3(odesl(:,1),odesl(:,2),odesl(:,3),'b')
plot3(soln(:,1),soln(:,2),soln(:,3),'r')
legend('ODE Solution','Circular Arc')
view([-30 -70])
hold off

Проблемы при установке ОДУ

Как описано в разделе Оптимизация моделирования или обычного дифференциального уравнения, оптимизатор может иметь проблемы из-за собственных шумов в числовых решениях ОДУ. Если вы подозреваете, что ваше решение не является идеальным, возможно, потому, что сообщение о выходе или флаг выхода указывают на потенциальную неточность, то попробуйте изменить конечное различие. В этом примере используйте больший размер шага конечной разности и центральные конечные разности.

options = optimoptions('lsqcurvefit','FiniteDifferenceStepSize',1e-4,...
    'FiniteDifferenceType','central');
[pbest2,presnorm2,presidual2,exitflag2,output2] = ...
    lsqcurvefit(@paramfun,xt0,tspan,soln,[],[],options);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

В этом случае использование этих вариантов конечных разностей не улучшает решение.

disp([presnorm,presnorm2])

   20.0637   20.0637