问题描述
我正在尝试绘制动态系统的状态空间图和时程图.不过有一个问题.状态空间被位于 x1 = 0 的平面分成两半.状态空间轴是 x1、x2、x3.x1 = 0 平面平行于 x2/x3 平面.x1 = 0 平面上方的状态空间由 eqx3 中的 ODE 描述，而 x1 = 0 平面下方的状态空间由 eqx4 中的 ODE 描述.
I am trying to plot a state-space diagram, as well as a time-history diagram, of a dynamical system. There's a catch, though. The state-space is divided into two halves by a plane located at x1 = 0. The state-space axes are x1, x2, x3. The x1 = 0 plane is parallel to the x2/x3 plane. The state-space above the x1 = 0 plane is described by the ODEs in eqx3, whereas the state-space below the x1 = 0 plane is described by the ODEs in eqx4.
所以，平面 x1 = 0 上存在不连续性.我有一个模糊的理解，即事件函数 (function [value,isterminal,direction] = myEventsFcn(t,y)) 应该可以使用，但我不知道给value"、isterminal"和direction"赋予什么值.
So, there is a discontinuity on the plane x1 = 0. I have a vague understanding that an event function (function [value,isterminal,direction] = myEventsFcn(t,y)) should be used, but I do not know what values to give to "value", "isterminal", and "direction".
在下面的代码中，我有一个 eqx3 的初始条件和另一个 eqx4 的初始条件.eqx3 的初始条件是在状态空间的上半部分 (x1 > 0).然后轨道撞到 x1 = 0 平面，出现不连续性，eqx4 轨迹从该平面开始，但与 eqx3 结束的点不同.
In my code below, I have an initial condition for eqx3, and another initial condition for eqx4. The initial condition for eqx3 is in the upper half of the state-space (x1 > 0). The orbit then hits the x1 = 0 plane and there is a discontinuity, and the eqx4 trajectory begins on this plane, but from a different point from where eqx3 ended.
这能做到吗?我如何将其放入代码中?当轨道到达平面 x1 = 0 时，我是否停止积分?
Can this be done? How do I put it into a code? Do I stop the integration when the orbit reaches the plane 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
推荐答案
无事件，使用近距离平滑系统
首先，由于系统之间的差异是常数向量的加法或减法，因此使用诸如 tanh 之类的 sigmoid 函数很容易找到系统的近似平滑版本.
Without events, using a close-by smooth system
First, as the difference between the systems is the addition or subtraction of a constant vector it is easy to find an approximate smooth version of the system using some sigmoid function like 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 和其他坐标接近零.
As is visible, after t=1.2 the solution is essentially constant with x(1)=0 and the other coordinates close to zero.
如果要使用事件机制，请使 ODE 和事件函数依赖于符号参数 S，表示过零的相位和方向.
If you want to use the event mechanism, make ODE and event function dependent on a sign parameter S denoting the phase and the direction of the zero crossing.
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.431.36 其中(理论上)x(1)=0 和矢量场从两侧指向另一个相位.理论上的解决方案是采用线性组合，将第一个分量设置为零，从而产生的方向与分离面相切.在上面的第一个变体中，sigmoid 自动实现了类似的功能.使用事件可以添加额外的事件函数来测试这些条件的向量场，以及当它们停止持续时.
The solution enters a sliding mode for 1.36 < t < 1.43 where (theoretically) x(1)=0 and the vector field points to the other phase from both sides. The theoretical solution is to take a linear combination that sets the first component to zero, so that the resulting direction is tangential to the separating surface. In the first variant above the sigmoid achieves something like this automatically. Using events one could add additional event functions that test the vector field for these conditions, and when they cease to persist.
一个快速的解决方案是加厚边界面，即测试x(1)+epsilon*S==0，这样在触发事件之前解必须穿过边界面.在滑动模式下，它会立即被推回，产生乒乓或之字形运动.epsilon 必须小以避免过多地干扰解决方案，但也不能太小以避免触发太多事件.使用 epsilon=2e-4 八度在滑动间隔的特写中给出以下解决方案
A quick solution is to thicken the boundary surface, that is, test for x(1)+epsilon*S==0 so that the solution has to cross the boundary surface before triggering the event. In sliding mode, it will be immediately pushed back, giving a ping-pong or zigzag motion. epsilon has to be small to not perturb the solution too much, however not too small to avoid the triggering of too many events. With epsilon=2e-4 octave gives the following solution in a close-up to the sliding interval
octave 求解器，在某种程度上也是 Matlab，如果它发生在第一个积分步骤中，将不会触发终端事件.出于这个原因，InitialStep 选项被设置为一个相当小的值，它应该约为 0.01*epsilon.完整的解决方案现在看起来类似于在第一个变体中获得的解决方案
The octave solver, and in some way also Matlab, will not trigger a terminal event if it happens in the first integration step. For that reason the InitialStep option was set to a rather small value, it should be about 0.01*epsilon. The full solution looks now similar to the one obtained in the first variant
这篇关于不确定如何在 Matlab 中使用事件函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！