2

我已经求解了一个初始值微分方程并绘制了 ODE45 给出的 Y 值。从图中我可以模糊地告诉根应该在哪里，但在给定的任务中，我需要非常准确地找到它。

我的第一个猜测是将多项式调整为我的 X 和 Y 值，然后求解多项式方程。但是我使用了 polyfit 并且有 69 个已知值，这给了我一个 68grade 的多项式，我无法解决。那么，有谁知道我如何在不知道实际方程的情况下找到一组给定 Y 值的“根”？任务中写到应该使用插值！

提前致谢！

4

2 回答 2

2

给定一个 Y 值向量（按照相应的 X 值稳步增加的意义进行排序），您可以很容易地找到根所在的 X 值附近。根要么在 Y 值为零的地方，要么在两个改变符号的连续 Y 值之间。此代码片段说明了这个想法：

X = -1:0.1:1;
Y = X.*X - 0.4;

root_exact_pos  = find(Y==0);
root_approx_pos = find(diff(sign(Y))~=0);

根在X值中，无论是 inX(root_exact_pos(k))还是介于X(root_approx_pos(k))and之间X(root_approx_pos(k)+1)，k从 1 到相应根位置数组的元素数。

从这里开始，您可以应用您想要找到更好的根近似值的任何插值（我会在 2 点之间使用线性）。

于 2013-04-13T13:11:34.277 回答
1

你说你需要“非常准确地”找到根。上面@CST-Link 的答案不是如果你用 数值求解微分方程怎么做ODE45。事实上，这是一个坏主意。它需要以高分辨率输出解点并引入错误，如果您正在积分的方程是保守的（即，您将有效地从总能量应始终保持恒定的解中添加或减去能量），这可能会特别糟糕。您需要使用 Matlab 的 ODE 求解器的事件检测（过零）功能。这通常能够以接近机器 epsilon 的精度找到您的根，eps. 而@CST-Link 给出的技术只会给出 ODE45 步长数量级的精度，这可能非常大（插值可用于帮助改善这一点，但eps除非你使用小步长）。

查看 Matlab 对 的帮助ODE45，事件功能似乎令人困惑，因此我将尝试使用基于ballode（help ballode以及edit ballode更多信息）的代码给出一个更简单的示例，Matlab 包含用于演示事件的示例。

function eventsdemo

options = odeset('Events',@efun); % specify name of events function
[t,y,te,ye,ie] = ode45(@f,[0 10],[0;20],options); % integrate
figure
plot(t,y(:,1),'b',te,ye(:,1),'r.')

function dydt = f(t,y)
dydt = [y(2);-9.8]; % differential equation for ballistic motion

function [value,isterminal,direction] = efun(t,y)
value = y(1);
isterminal = 1;
direction = -1;

这个例子只是一个球在垂直方向上的简单弹道运动：函数f。目标是准确检测球在穿过地平面时的速度，y(1) == 0线在负方向。如果你过早发现接触，球的速度将低于实际速度，能量将丢失；为时已晚，能量将被注入。的value输出efun定义了交叉点必须为零的方程。这可以根据需要简单或复杂，可以取决于所有状态变量和时间，并且您可以通过指定向量来检测多个交叉点。在您的情况下，听起来您对 x 轴上的根感兴趣，所以我想您可能有相同的定义value. 如果您只想要第一个根或者只有一个，则isterminal在找到它时将停止集成。最后，如果您不知道函数在根部的斜率/梯度，您可以设置direction为零。用上面的代码试试这个事件函数（ie可以用来判断在te和ye输出中触发了三个事件中的哪一个）：

function [value,isterminal,direction] = efun(t,y)
value = [y(2) y(1)-10 y(1)];
isterminal = [0 0 1];
direction = [-1 0 -1];

一些小警告。您需要将最终积分时间设置为足够长。换句话说，事件必须在t0and之间发生（如果您不知道事件将在哪里发生，tf请参阅ballode迭代调用的方案）。ODE45如果设置isterminal为零，则输出t向量和y矩阵将被修剪以在终端事件发生时结束。最后，如果您指定固定步长输出，例如 ，TSPAN = t0:dt:tf最后一个输出点的步长可能会更小。

于 2013-04-13T19:51:28.273 回答