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

Question

我正在尝试绘制动态系统的状态空间图以及时间历史图。不过，有一个问题。状态空间被位于 x1 = 0 的平面分成两半。状态空间轴是 x1、x2、x3。 x1 = 0 平面平行于 x2/x3 平面。 x1 = 0 平面上方的状态空间由 eqx3 中的 ODE 描述，而 x1 = 0 平面下方的状态空间由 eqx4 中的 ODE 描述。

所以，平面 x1 = 0 上存在不连续性。我有一个模糊的理解，应该使用事件函数（function [value,isterminal,direction] = myEventsFcn(t,y)），但我不知道要给“value”、“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 之类的 sigmoid 函数很容易找到系统的近似平滑版本。

  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 并且矢量场从两侧指向另一个相位。理论上的解决方案是采用线性组合，将第一个分量设置为零，这样得到的方向与分离面相切。在上面的第一个变体中，sigmoid 自动实现了类似的功能。使用事件可以添加额外的事件函数来测试矢量场的这些条件，以及它们何时停止存在。

一个快速的解决方案是加厚边界表面，即测试x(1)+epsilon*S==0，以便解决方案必须在触发事件之前穿过边界表面。在滑动模式下，它会立即被推回，做出乒乓或曲折的运动。 epsilon 必须很小，以免过多干扰解决方案，但也不能太小，以免触发太多事件。用epsilon=2e-4 octave 在滑动区间的特写中给出了以下解

如果发生在第一个积分步骤中，倍频程求解器，以及在某种程度上还包括 Matlab，将不会触发终端事件。由于这个原因，InitialStep 选项被设置为一个相当小的值，它应该是大约0.01*epsilon。完整的解决方案现在看起来类似于在第一个变体中获得的解决方案