不确定如何在Matlab中使用事件函数

Question

我正在尝试绘制动力系统的状态空间图以及时间历史图。不过有一个陷阱。状态空间被位于x1 = 0的平面分为两半。状态空间轴为x1，x2，x3。 x1 = 0平面平行于x2 / x3平面。 x1 = 0平面上方的状态空间由eqx3中的ODE描述，而x1 = 0平面下方的状态空间由eqx4中的ODE描述。因此，平面x1 = 0上存在一个不连续点。我模糊地理解应该使用事件函数（函数[value，isterminal，direction] = myEventsFcn（t，y）），但是我不知道什么值赋予“值”，“ isterminal”和“方向”。在下面的代码中，我对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之类的S形函数来找到系统的近似平滑版本。

  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);
 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)myEventsFcn1(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] = myEventsFcn1(t,x,S)
   value = x(1);
   isterminal = 1;
   direction = -S;
 end
 function [value,isterminal,direction] = myEventsFcn2(t,x,S)
   value = x(1)+S*1e-8;
   isterminal = 1;
   direction = 0;
 end
我无法使用GNU八度获得一致的结果，因为未发现x(1)的一些明显的零交叉。作为一种解决方法，我还尝试了第二个事件函数，该事件在进入相反符号的区域后不久但确定地触发了该事件。但是结果又令人难以置信。使用Matlab可能会获得更好的成功。第二个变体中的八度给出了图