ODE with time dependent input, How to speed Up without using interpolation?

Question

I am trying to solve a system of ODEs and my input excitation is a function of time.

I have been using interp1 inside the integration function, but this doesn't seems like a very efficient way to do this. I know it is not, because once I change the input excitation to a sin function, which does not require an interp1 call inside the function, I get much much faster results. But doing interpolation every step takes about 10–20 times longer to converge. So, is there a better way of solving ODEs for arbitrary time dependent excitation, without needing to do interpolation or some other tricks to speed up?

I am just copying a modified version of a simple example from The MathWorks here:

Input Excitation is a gradually increasing sin function, but after some time later it becomes a constant amplitude sin function.

Dt   = 0.01;   % sampling time step
Amp0 = 2;      % Final Amplitude of signal
Dur_G = 10;    % Duration of gradually increasing part of signal 
Dur_tot = 25;  % Duration of total signal
t_G = 0 : Dt : Dur_G; % time of gradual part 
A = linspace(0, Amp0, length(t_G));
carrier_1 = sin(5*t_G); % Unit Normal Signal
carrier_A0 = Amp0*sin(5*t_G);
out_G = A.*carrier_1;   % Gradually Increasing Signal
% Total Signal with Gradual Constant Amplitude Parts
t_C = Dur_G+Dt:Dt:Dur_tot; % time of constant part
out_C = Amp0*sin(5*t_C);   % Signal of constant part
ft = [t_G t_C];    % total time 
f = [out_G out_C]; % total signal
figure; plot(ft, f, '-b'); % input excitation


function dydt = myode(t,y,ft,f)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = 2;               % a constant
dydt = -f.*y + g;    % Evaluate ODE at time t

tspan = [1 5]; ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f), tspan, ic, opts);

figure;
plot(t,y);

Note that I explained only first part of my problem above, which is solving system for a gradually increasing sin function. In the second part, I need to solve it for an arbitrary input excitation (e.g., a ground acceleration input).

Answer 1

For this example, you could use griddedInterpolant class to get a bit of a speed-up:

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;
Fp = griddedInterpolant(ft,f);
gt = linspace(1,6,25);
g = 3*sin(gt-0.25);
Gp = griddedInterpolant(gt,g);

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y)myode(t,y,Fp,Gp),tspan,ic,opts);

figure;
plot(t,y);

The ODE function is then:

function dydt = myode(t,y,Fp,Gp)
f = Fp(t);        % Interpolate the data set (ft,f) at time t
g = Gp(t);        % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t

On my system with R2015b, the call to ode45 is about three times faster (0.011 sec vs. 0.035 sec) for your example. You could get a bit more speed by switching to ode23. You can read more about the griddedInterpolant class here.

If your actual system, discretely switches between inputs particular points in time, then you should probably solve the problem piecewise by integrating each case separately. See this question and this question. If the system switches based on the value of the state variable(s), then you should use event location (see this question). However, if "solving ODEs for random time dependent excitation" means that you're adding random noise to the system, then you have an SDE rather than an ODE, which is a completely different beast.