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

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

Значения Лоренца параметров для чувствительной системы $$\sigma =10,\beta =8/3,\rho =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

Соответствуйте ОДУ к круговой дуге

Теперь измените параметры $$\sigma ,\beta ,and\rho$$ лучше всего соответствовать круговой дуге. Для еще лучшей подгонки позвольте начальной точке [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