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

Question

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

所以，x1 = 0平面上有一个不连续点。我有一个模糊的理解，应该使用事件函数（function [value,isterminal,direction] = myEventsFcn(t,y)），但我不知道给“值”赋什么值"、"isterminal" 和 "direction"。

在我下面的代码中，我有一个 eqx3 的初始条件，以及另一个 eqx4 的初始条件。 eqx3 的初始条件在状态的上半部分-space (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

没有事件，使用附近的平滑系统

首先，由于系统之间的区别在于常数向量的加法或减法，因此使用一些 S 形函数很容易找到系统的近似平滑版本，例如 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 并且向量场从两侧指向另一相。理论上的解决方案是采用将第一个分量设置为零的线性组合，以便生成的方向与分离表面相切。在上面的第一个变体中，sigmoid 自动实现了类似的功能。使用事件可以添加额外的事件函数来测试矢量场的这些条件，以及它们何时停止持续存在。

一个快速的解决方案是加厚边界表面，即测试x(1)+epsilon*S==0，使得解决方案必须穿过边界表面才能触发事件。在滑动模式下，它会立即被推回，进行乒乓或之字形运动。 epsilon 必须很小以免对解决方案造成太大干扰，但也不能太小以避免触发太多事件。使用 epsilon=2e-4 octave 在滑动间隔

的特写中给出以下解决方案

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

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

Unsure about how to use event function in Matlab

matlab

numerical-methods

ode

没有事件，使用附近的平滑系统

有活动