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

Система Лоренца является системой обыкновенных дифференциальных уравнений (см. Система Лоренца). Для действительных констант $$\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