插值后用ode45求解ODE

Question

我正在处理以下代码，但找不到问题所在：

xx = 0:1/50:1;
v = 3.*exp(-xx)-0.4*xx;
xq = xx;
vq = @(xq) interp1(xx,v,xq);

tspan = 0:1/50:1;
x0 = 3;
[~, y2] = ode45(@(t,x)vq(x), tspan, x0);

我明白y2 = [3;NAN;NAN;NAN,.....]的意思。然而，当我在调用ode45之前绘制这两个方程时，我会发现它们是相等的，这并不奇怪。

当我计算：

f = @(t,r) 3.*exp(-r)-0.4*r;
[~, y] = ode45(f,tspan,x0);

效果很好。但是我需要证明，如果我插值，我可以得到同样的结果。为什么不起作用？

Answer 1

您可以获得NaN，因为这是interp1为xx跨越的间隔以外的值返回的默认值。在这种情况下，xx只会因0而异于1。但你的初始条件是在3。如果要使用插值，则需要从数据定义的间隔内开始，并确保保持在该区间内。例如，如果您只是更改初始条件：

xx = 0:1/50:1;
v = 3.*exp(-xx)-0.4*xx;
xq = xx;
vq = @(xq) interp1(xx,v,xq);

tspan = 0:1/50:1;
x0 = 0.1;
[t, y2] = ode45(@(t,x)vq(x), tspan, x0);

f = @(t,r) 3.*exp(-r)-0.4*r;
[t, y] = ode45(f,tspan,x0);

figure;
subplot(211)
plot(t,y,'b',t,y2,'r--')
subplot(212)
plot(t,abs(y-y2))
xlabel('t')
ylabel('Absolute Error')

即使在这个初始条件下，由于指数增长，在某个点上，系统的状态留下您的0，1间隔，y2将成为NaN。不过，如果您愿意的话，可以告诉interp1使用实际外推。