How can I make ODE15s interpolate a value from a given table on each time step?

I am trying to solve a system of algebraic and differntial equations using ODE15s. This is not my problem though, I cannot seem to figure out how to impliment my experimental data, which I have in a matrix in terms of time, into my ODE solver. I would basically like for, on each time step, my Area of Valve(Av) in the ODE solver to interpolate to a matrix containing Valve area and time that I have measured.

So just for an example here's a much simplified version of my code:

Data = readmatrix(experimentdata)
Av = Data(:,8) %Table of Av values that I should note are too sporatic to be polyfit with any accuracy.
time = Data(:,1)

function out = myODE(t,y)
out(1) = ...
out(2) = ...
out(3) = Av*...
end

Basically because I can't get my Av value in terms of time to match with the time step in the ODE solver, I am getting tolerance errors from ODE15s and having discontinuities created. Any help is greatly appreciated.

Note: The code runs just fine given constant Av, only when given a vector does it fail.

Answer 1

[Walter Roberson]

You can call interp1(), but you must use an interpolation method that is at least quartic (not linear).

If your situation is such that the arrays you are interpolating over can be said to outline a continuous path that should be assumed to curve through the data points, then interp1() with 'cubic' might work for you. If you are at time t between two entries t1 and t2 and it makes sense to "look ahead" to t2 to steer towards where you need to be at t2, then this method might work for you.

If your situation is such that the arrays give information about sudden events, then each time there is such an event, then you need to stop the ode*() call, adjust the boundary conditions, and then call the ode*() function again to pick up from the boundary. For example if the data in the arrays represent injections of chemicals, or represent bumps in a road, or represent bouncing off of something, then you cannot use interp1() in your ode function.

Answer 2

[Star Strider]

That depends on the problem that you want ode15s to solve. As a general rule, the approach is to include an interp1 call, and pass ‘Av’ and ‘time’ as additional arguments to ‘myODE’.

Example —

time = linspace(0, 5, 150).';
Av = (1 - exp(-0.5*time))*10;

[t,y] = ode15s(@(t,y)myODE(t,y,Av,time), [0 5], [0.1 0.2 0.3]);

figure
plot(t, y)
grid
xlabel('Time')
ylabel('Something')

function out = myODE(t, y, Av, time)
Avv = interp1(Av, time, t);
out = zeros(3,1);
out(1) = -100*y(1);
out(2) =  y(2);
out(3) = Avv / y(3);
end

That may not be appropriate for all types of problems that ode15s can solve, however it will work for some of them.

.

[Walter Roberson]

Never use the default linear interpolation of interp1() within any of the ode* functions; doing so violates the mathematical requirements of the Runge-Kutta calculations.

[Star Strider]

I used the equivalent of the ode45 documentation example in ODE with Time-Dependent Terms. Other than fitting a nonlinear regression model to the ‘Av’ data and then creating a predict function from it, I see no other option.

[Walter Roberson]

If you use the default linear interpolation for interp1() inside the ode function, then you must either use event functions or else you must make careful use of the time boundaries, so that there is no call to the ode*() routines that crosses one of the times for which there is data. (Though you could get away with it if the function being approximated is constant or linear.)

[Torsten]

@Walter Roberson

Never use the default linear interpolation of interp1() within any of the ode* functions; doing so violates the mathematical requirements of the Runge-Kutta calculations.

In theory, the functions that define the ODEs have to be differentiable as often as the order of the integration method used. But in practice, this is rarely true.

[Noah Cain]

I have attempted this solution without success. Essentially adding an interp1 call before the differential equations as above, however this time I used spline since previously I had been using regular interp1 and would receive error messages about integration tolerances. However even in using spline as below I still am receiving the same error(also below)

Av = interp1(Valvekept,Timekept,t,'spline') ;

Warning: Failure at t=0.000000e+00.  Unable to meet integration
tolerances without reducing the step size below the smallest value
allowed (7.905050e-323) at time t. 
> In ode15s (line 661)
In ode15_20D_withPID (line 51) 

If it helps a all the data in Valvekept which i'm trying to interpolate from looks like this:

[Star Strider]

The integration tolerance Warning may have nothing to do with the interpolation, although it may have something to do with the values being interpolated, that being the results of the interpolation rather than the interpolation itself. (Note that extrapolating, for example if the time being interpolated is not with in the range of the values of time corresponding to the ‘Av’ values, and the function is not told to extrapolate, the result of the interpolation will be NaN. This will likely not produce a warning, however it will produce NaN values for that differential equation and all others that depend on it, for the rest of the integration.)

The Warning means that the system of differential equations is producing a value of ±Inf, and that causes the integration to stop. Since I do not know any of the details, I cannot determine what the problem is.

[Torsten]

@Noah Cain

The call to "interp1" should most probably be

Av = interp1(Timekept,Valvekept,t,'spline') ;

instead of

Av = interp1(Valvekept,Timekept,t,'spline') ;

shouldn't it ?

[Star Strider]

@Torsten — Thank you. I didn’t see that (but then I wasn’t looking for it).

That doesn’t affect my previous Comment, that remains the same.