事件发生时Matlab ode45不终止
Matlab ode45 Does Not Terminate When Event Occurs
当使用 MATLAB 的 ode45
求解器求解 ODE 时,我试图触发一个事件以在满足(一般)布尔条件时终止求解。例如,在下面的代码中,条件是 x > 1.5
(因为事件在值为零时触发,即当 x <= 1.5
为假时)。从条件在时间零时满足但之后没有立即满足的点开始解决方案会导致永远不会触发事件。这是怎么回事?这是过零检测中的错误吗?我知道使用连续值时过零更容易,但我希望求解器在这种情况下仍然可以工作,尽管速度较慢。
f = @(t, x) 1;
tspan = [0, 10];
x0 = 1.5;
[t, x] = ode45(f, tspan, x0, odeset('Events',@zeroevents));
x(end) % Value is 11.5
assert(x(end) <= 1.5001) % This assertion fails.
function [value,isterminal,direction] = zeroevents(t, x)
isterminal = 1;
direction = 0;
value = x <= 1.5;
end
此行为可能是由于以下三个原因中的两个 limitations of event detection:
- If a terminal event occurs during the first step of the integration, then the solver registers the event as nonterminal and continues integrating.
和
- Zeros are determined by sign crossings between steps. Therefore, zeros of functions with an even number of crossings between steps can be missed.
因此,尽管初始条件满足要求,但根据定义忽略了第一次触发事件。由于您的值是通过布尔表达式计算的,因此它永远不会出现符号变化,因此永远不会触发检测。
将您的问题更改为在初始步骤后触发事件并通过带符号的残差(即 value = x - 1.5;
)计算值应该会产生您正在寻找的行为。如果接近值的方式很重要(即从 above/positive 或 below/negative 开始),则应适当设置方向值;两种方法的零火方向。
我的解决方法是在解决方案中第二次手动检查是否满足停止条件。如果不是,那么我将截断解决方案的其余部分。
f = @(t, x) 1;
tspan = [0, 10];
x0 = 1.5;
[t, x] = ode45(f, tspan, x0, odeset('Events',@zeroevents));
x(end) % Value is 11.5
% Truncate the rest of the solution.
if isStop(x(2))
t = t(1:2);
x = x(1:2);
end
x(end) % Value is ~1.57
function [value,isterminal,direction] = zeroevents(t, x)
isterminal = 1;
direction = 0;
value = isStop(x);
end
function stop = isStop(x)
stop = x > 1.5;
end
然而,这不是一个很好的解决方案,因为 ode45
选择了一个相对较大的时间步长,因为它没有“注意到”事件,所以 x(2)
明显大于 1.5
(大约 1.57
)。要解决此问题,我们可以使用具有小时间步长的线性外推法手动计算 x(2)
:
if isStop(x(2))
delta_t = 1e-10;
t1 = t(1);
t2 = t1 + delta_t;
t = [t1; t2];
x1 = x(1);
x2 = x1 + delta_t * f(x1, t1);
x = [x1'; x2'];
end
当使用 MATLAB 的 ode45
求解器求解 ODE 时,我试图触发一个事件以在满足(一般)布尔条件时终止求解。例如,在下面的代码中,条件是 x > 1.5
(因为事件在值为零时触发,即当 x <= 1.5
为假时)。从条件在时间零时满足但之后没有立即满足的点开始解决方案会导致永远不会触发事件。这是怎么回事?这是过零检测中的错误吗?我知道使用连续值时过零更容易,但我希望求解器在这种情况下仍然可以工作,尽管速度较慢。
f = @(t, x) 1;
tspan = [0, 10];
x0 = 1.5;
[t, x] = ode45(f, tspan, x0, odeset('Events',@zeroevents));
x(end) % Value is 11.5
assert(x(end) <= 1.5001) % This assertion fails.
function [value,isterminal,direction] = zeroevents(t, x)
isterminal = 1;
direction = 0;
value = x <= 1.5;
end
此行为可能是由于以下三个原因中的两个 limitations of event detection:
- If a terminal event occurs during the first step of the integration, then the solver registers the event as nonterminal and continues integrating.
和
- Zeros are determined by sign crossings between steps. Therefore, zeros of functions with an even number of crossings between steps can be missed.
因此,尽管初始条件满足要求,但根据定义忽略了第一次触发事件。由于您的值是通过布尔表达式计算的,因此它永远不会出现符号变化,因此永远不会触发检测。
将您的问题更改为在初始步骤后触发事件并通过带符号的残差(即 value = x - 1.5;
)计算值应该会产生您正在寻找的行为。如果接近值的方式很重要(即从 above/positive 或 below/negative 开始),则应适当设置方向值;两种方法的零火方向。
我的解决方法是在解决方案中第二次手动检查是否满足停止条件。如果不是,那么我将截断解决方案的其余部分。 f = @(t, x) 1; tspan = [0, 10]; x0 = 1.5; [t, x] = ode45(f, tspan, x0, odeset('Events',@zeroevents));
x(end) % Value is 11.5
% Truncate the rest of the solution.
if isStop(x(2))
t = t(1:2);
x = x(1:2);
end
x(end) % Value is ~1.57
function [value,isterminal,direction] = zeroevents(t, x)
isterminal = 1;
direction = 0;
value = isStop(x);
end
function stop = isStop(x)
stop = x > 1.5;
end
然而,这不是一个很好的解决方案,因为 ode45
选择了一个相对较大的时间步长,因为它没有“注意到”事件,所以 x(2)
明显大于 1.5
(大约 1.57
)。要解决此问题,我们可以使用具有小时间步长的线性外推法手动计算 x(2)
:
if isStop(x(2))
delta_t = 1e-10;
t1 = t(1);
t2 = t1 + delta_t;
t = [t1; t2];
x1 = x(1);
x2 = x1 + delta_t * f(x1, t1);
x = [x1'; x2'];
end