Не знаете, как использовать функцию события в Matlab

Question

Я пытаюсь построить диаграмму пространства состояний, а также диаграмму истории времени динамической системы. Хотя есть подвох. Пространство состояний разделено на две половины плоскостью, расположенной при x1 = 0. Оси пространства состояний имеют вид x1, x2, x3. Плоскость x1 = 0 параллельна плоскости x2 / x3. Пространство состояний над плоскостью x1 = 0 описывается ODE в eqx3, тогда как пространство состояний ниже плоскости x1 = 0 описывается ODE в eqx4. Таким образом, на плоскости x1 = 0 есть разрыв дать «значение», «истерминал» и «направление». В моем коде ниже у меня есть начальное условие для eqx3 и другое начальное условие для eqx4. Начальное условие для eqx3 находится в верхней половине пространства состояний (x1> 0). Затем орбита попадает в плоскость x1 = 0, и возникает разрыв, и траектория eqx4 начинается на этой плоскости, но не в той точке, где закончился eqx3. Можно ли это сделать? Как мне вставить это в код? Останавливаю ли я интегрирование, когда орбита достигает плоскости x1 = 0?

eta = 0.05;
omega = 25;
tspan = [0,50];
initcond = [2, 3, 4]
[t,x] = ode45(@(t,x) eqx3(t,x,eta, omega), tspan, initcond);
initcond1 = [0, 1, 1]
[t,y] = ode45(@(t,y) eqx4(t,y,eta, omega), tspan, initcond1);

plot3(x(:,1), x(:,2), x(:,3),y(:,1), y(:,2), y(:,3))
xlabel('x1')
ylabel('x2')
zlabel('x3')

%subplot(222)
%plot(t, x(:,1), t,x(:,2),t,x(:,3),'--');
%xlabel('t')

function xdot = eqx3(t,x,eta,omega)
  xdot = zeros(3,1);
  xdot(1) = -(2*eta*omega + 1)*x(1) + x(2) - 1;
  xdot(2) = -(2*eta*omega + (omega^2))*x(1) + x(3) + 2;
  xdot(3) = -(omega^2)*x(1) + x(2) - 1;

end

function ydot = eqx4(t,y,eta,omega)
  ydot = zeros(3,1);
  ydot(1) = -(2*eta*omega + 1)*y(1) + y(2) + 1;
  ydot(2) = -(2*eta*omega + (omega^2))*y(1) + y(3) - 2;
  ydot(3) = -(omega^2)*y(1) + y(2) + 1;

end

function [value,isterminal,direction] = myEventsFcn(t,y)
   value = 0
   isterminal = 1
   direction = 1

end

Answer 1

Без событий, используя близкую гладкую систему

Во-первых, поскольку различие между системами заключается в сложении или вычитании постоянного вектора, легко найти приближенную гладкую версию системы, используя некоторые сигмовидная функция, такая как tanh.

  eta = 0.05;
  omega = 25;
  t0=0; tf=4;
  initcond = [2, 3, 4];
  opt = odeset('AbsTol',1e-11,'RelTol',1e-13);
  opt = odeset(opt,'MaxStep',0.005); % in Matlab: opt = odeset('Refine',10);
  [t,x] = ode45(@(t,x) eqx34(t,x,eta, omega), [t0, tf], initcond, opt);

  clf;
  subplot(121);
  plot3(x(:,1), x(:,2), x(:,3));
  xlabel('x1');
  ylabel('x2');
  zlabel('x3');

  subplot(122);
  plot(t, 10.*x(:,1), t,x(:,2),':',t,x(:,3),'--');
  xlabel('t');

  function xdot = eqx34(t,x,eta,omega)
    S = tanh(1e6*x(1));
    xdot = zeros(3,1);
    xdot(1) = -(2*eta*omega + 1)*x(1) + x(2) - S;
    xdot(2) = -(2*eta*omega + (omega^2))*x(1) + x(3) + 2*S;
    xdot(3) = -(omega^2)*x(1) + x(2) - S;
  end

, в результате чего графики

Как видно, после t=1.2 решение по существу постоянно с x(1)=0, а другие координаты близки к нулю.

С событиями

Если вы хотите использовать механизм событий, сделайте ODE и функцию события зависимыми от параметра знака S обозначает фазу и направление пересечения нуля.

 eta = 0.05;
 omega = 25;
 t0=0; tf=4;
 initcond = [2, 3, 4];
 opt = odeset('AbsTol',1e-10,'RelTol',1e-12,'InitialStep',1e-6);
 opt = odeset(opt,'MaxStep',0.005); % in Matlab: opt = odeset('Refine',10);
 T = [t0]; TE = [];
 X = [initcond]; XE = [];
 S = 1; % sign of x(1)
 while t0<tf
    opt = odeset(opt,'Events', @(t,x)myEventsFcn(t,x,S));
    [t,x,te,xe,ie] = ode45(@(t,x) eqx34(t,x,eta, omega, S), [t0, tf], initcond, opt);
    T = [ T; t(2:end) ]; TE = [ TE; te ];
    X = [ X; x(2:end,:)]; XE = [ XE; xe ];
    t0 = t(end);
    initcond = x(end,:);
    S = -S % alternatively = 1-2*(initcond(1)<0);
 end
 disp(TE); disp(XE);

 subplot(121);
 hold on;
 plot3(X(:,1), X(:,2), X(:,3),'b-');
 plot3(XE(:,1), XE(:,2), XE(:,3),'or');
 hold off;
 xlabel('x1');
 ylabel('x2');
 zlabel('x3');

 subplot(122);
 plot(T, 10.*X(:,1), T,X(:,2),':',T,X(:,3),'--');
 xlabel('t');

 function xdot = eqx34(t,x,eta,omega,S)
  xdot = zeros(3,1);
  xdot(1) = -(2*eta*omega + 1)*x(1) + x(2) - S;
  xdot(2) = -(2*eta*omega + (omega^2))*x(1) + x(3) + 2*S;
  xdot(3) = -(omega^2)*x(1) + x(2) - S;
 end

 function [value,isterminal,direction] = myEventsFcn(t,x,S)
   value = x(1)+2e-4*S;
   isterminal = 1;
   direction = -S;
 end

Решение входит в скользящий режим для 1.36 < t < 1.43, где (теоретически) x(1)=0, а векторное поле указывает на другую фазу из обоих стороны. Теоретическое решение состоит в том, чтобы взять линейную комбинацию, которая устанавливает первый компонент на ноль, так что результирующее направление является касательным к разделяющей поверхности. В первом варианте выше сигмоид достигает чего-то подобного автоматически. Используя события, можно добавить дополнительные функции событий, которые проверяют векторное поле на эти условия и когда они перестают существовать.

Быстрое решение состоит в том, чтобы утолщить граничную поверхность, т. Е. Проверить на x(1)+epsilon*S==0, так что решение должно пересечь граничную поверхность, прежде чем вызвать событие. В скользящем режиме он будет сразу же сдвинут назад, давая пинг-понг или зигзагообразное движение. epsilon должно быть маленьким, чтобы не слишком мешать решению, но не слишком маленьким, чтобы избежать запуска слишком большого количества событий. С epsilon=2e-4 Октава дает следующее решение в приближении к скользящему интервалу

Октава решатель, а в некотором роде и Matlab не будет запускать терминальное событие, если оно происходит на первом этапе интеграции. По этой причине для параметра InitialStep было установлено довольно небольшое значение, оно должно быть около 0.01*epsilon. Полное решение теперь похоже на решение, полученное в первом варианте