Соответствуйте обыкновенному дифференциальному уравнению (ODE)

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

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

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

$\begin{array}{l} \frac{d x}{d t} = σ (y - x) \\ \frac{d y}{d t} = x (ρ - z) - y \\ \frac{d z}{d t} = x y - β z . \end{array}$

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

Соответствуйте круговому пути к решению для ОДУ

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

\angle $θ (1)$ из пути от x-y плоскости
\angle $θ (2)$ из плоскости от наклона вдоль оси X
Радиус R
Скорость V
Переключите t0 со времени 0
3-D сдвиг в дельте пробела

С точки зрения этих параметров определите положение кругового пути в течение многих времен 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

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

Теперь измените параметры $σ, β, a n d ρ$ лучше всего соответствовать круговой дуге. Для еще лучшей подгонки позвольте начальной точке [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

Документация