How to pass step forcing function in system model to ode15s

I am using ode15s to solve a system that has quasi-step forcing inputs that are defined by data. I have used griddedinterpolant with a variety of methods but all result in a smoothed forcing that is not realistic for the problem. The data are available in time steps of 10 minutes, and I know that at minute 30, the infusions change abruptly. How can I represent the forcings to ode15s to account for the fact that it has a variable time step and needs to compute the intermediate data?

Answer 1

[Torsten]

You will have to stop the integration when the inputs get discontinuous and restart with the results obtained so far, maybe with new initial conditions and new equations.

Answer 2

[Sam Chak]

Hi @Andrew

While the ODE solver can change its integration step size on the fly to find sharp corners, by drastically shrinking its step size to , the gold standard for handling discontinuities in numerical simulation is to stop right at the sharp corner, reset, and restart with a clean slate.

If you know exactly when the discontinuities happen, you can use the stopping-and-restarting method (in just 5 steps) so that the solver integrates smoothly up to the exact time step of the change, stops, and starts fresh without carrying over historical error terms.

Here is the basic code setup:

%% STEP 1: Visually inspect the External input signal
tau     = 20;       % period of signal
Tf      = 40;       % duration of signal
Ts      = 0.01;     % step size of signal
[u, tu] = gensig("square", tau, Tf, Ts);
figure
plot(tu, u), grid on
ylim([-1, 2])
title('External input signal u(t) with discontinuous jumps')
xlabel('Time (minutes)')
ylabel('Amplitude')

%% STEP 2: Identify the jumping segments where the discontinuities happen 
T   = 0:10:40

T = 1×5

0 10 20 30 40

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% STEP 3: Test the loop first
figure
hold on
for k = 1:(length(T) - 1)
    tspan   = [T(k), T(k+1)]
    idx     = tu >= T(k) & tu <= T(k+1);
    plot(tu(idx), u(idx));
end

tspan = 1×2

0 10

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

tspan = 1×2

10 20

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

tspan = 1×2

20 30

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

tspan = 1×2

30 40

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

hold off
grid on
ylim([-1, 2])
title('Plot external input signal u(t) after each jumping segment')
xlabel('Time (minutes)')
ylabel('Amplitude')

%% STEP 4: Describe the ODE with time-dependent parameters
function dx = system45(t, x, tu, u)
    % External input forced signal, u(t)
    u       = interp1(tu, u, t);    % Interpolate the data set (tu, u) at time t
    
    % Dynamics of system45
    dx      = zeros(2, 1);
    dx(1)   = x(2);
    dx(2)   = 4*u - 4*x(1) - 4*x(2);
end

%% STEP 5: Loop the ODE calls and plot the result simultaneously
x0  = [1; 0];   % Initial conditions at t=0

figure
hold on
for k = 1:(length(T) - 1)
    % identify simulation segment before the next jump
    tspan   = [T(k), T(k+1)];
    idx     = tu >= T(k) & tu <= T(k+1);
    
    % run ODE solver (for each tspan segment)
    [t, x]  = ode45(@(t, x) system45(t, x, tu(idx), u(idx)), tspan, x0);
    
    % plot the 1st component of the solution
    plot(t, x(:, 1))
    
    % extract final state to use as the NEW initial condition
    x0      = x(end, :);
end
hold off
grid on
ylim([-1, 2])
title('Time history of the 1st component of the solution')
xlabel('Time (minutes)')
ylabel('Amplitude')