1

我有以下简单的 ODE：

dx/dt=-1

初始条件 x(0)=5，我对 x(t)==1 感兴趣。所以我有以下事件功能：

function [value,isterminal,direction] = test_events(t,x)
    value = x-1;
    isterminal = 0;
    direction = 0;
end

这应该在 t=4 时产生一个事件。但是，如果我运行以下代码，我会得到两个事件，一个在 t=4，一个在附近的位置 t=4+5.7e-14：

options = odeset('Events',@test_events);
sol = ode45(@(t,x)-1,[0 10],5,options);
fprintf('%.16f\n',sol.xe)
% 4.0000000000000000
% 4.0000000000000568

如果我运行类似的代码来查找 x(t)==0 或 x(t)==-1 （分别为 value = x; 或 value = x+1; ），我只有一个事件。为什么这会产生两个事件？

更新：如果选项结构更改为以下内容：

options = odeset('Events',@test_events,'RelTol',1e-4);

...那么 ODE 仅在 t=4+5.7e-14 返回一个事件。如果 'RelTol' 设置为 1e-5，它会在 t=4 时返回一个事件。如果 'RelTol' 设置为 1e-8，它返回与默认相同的两个事件 ('RelTol'=1e-3)。此外，将初始条件从 x(0)=5 更改为 x(0)=4 会产生一个事件，但设置 x(0)=4 和 'RelTol'=1e-8 会产生两个事件。

更新 2：观察 sol.x 和 sol.y 输出（分别为 t 和 x），时间以整数 [0 1 2 3 4 5 6 7...] 进行，并且 x 作为整数进行直到 x(t =5) 像这样：[5 4 3 2 1 1.11e-16 -1.000 -2.000...]。这表明在 t=4 和 t=5 之间发生了一些事情，在 ODE 解决方案中产生了“凹凸”。为什么？

4

1 回答 1

0

一种推测可以解释在这个简单问题中如何发生舍入误差：使用k_nODE 导数函数的评估在内部步骤之间插入解决方案，也称为“密集输出”。理论形式为

b_1(u)k_1 + b_2(u)k_2 + ...b_s(u)k_s

其中0 <= u<= 1它是内部点之间间隔上的参数，即t = (1-u)*t_k+u*t_{k+1}.

系数多项式是非平凡的。虽然在示例中所有k_i=1都是常数，但求和的计算b_1(u)+...+b_s(u)会累积舍入误差，这些舍入误差在接近根的解值中变得可见，即使y_k和y_{k+1}是精确的。在该累积浮点噪声范围内，该值可能会围绕根部振荡，从而导致检测到多个过零。

于 2019-03-07T16:55:08.343 回答