Solving an ODE with ode45 after interpolating

Question

I am working with the following code and can't find what is wrong:

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);

I get that y2 = [3;NAN;NAN;NAN,.....]. Yet, when I plot both equations before calling ode45, I get that they are equal, which is not a surprise.

When I compute:

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

it works fine. But I need to show that I can get the same results if I interpolate. Why is that not working?

Answer 1

You get NaN because that is the default returned by interp1 for values outside of the interval spanned by xx. In you case, xx only varies from 0 to 1. But your initial condition is at 3. If you're going to use interpolation you need start inside of the interval defined by your data and make sure you stay there. For example, if you simply change your initial condition:

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')

Even with this initial condition, because of the exponential growth, at a certain point the state of the system leaves your [0, 1] interval and y2 will become NaN. You can tell interp1 to use actual extrapolation if you prefer though.