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.

[Andrew]

Thanks--this may be difficult to do when the past history precedes the breakpoint.

[Steven Lord]

Can you say a little more about the differential equation you're trying to solve? Is it possible you're not trying to solve an ordinary differential equation (for which ode15s might be the right tool) but you're trying to solve a delay diferential equation (for which you may want dde23 or one of its "cousins")?

[Andrew]

Hi Steven,

I'm modeling a 2-compartment mass transfer system described by 4 state variables and forced by one or two external imputs depending upon the experimental situation modeled. It has been traditionally modeled as a set of ODEs and others and I have had success using ODE15s. I'm unfamiliar with the applicaiton of delay differential approaches but will look them up.

Thanks for your evaluation and advice.

Andy

[Sam Chak]

Hi Andy,

For your info, a delayed input in a dynamic system is fundamentally different from a delay-differential equation (DDE).

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

[Andrew]

Thanks for your detailed eaxamination and explanation! I think I'll try resampling the forcing functions to a shorter dt and use that shorter dt in the system solution.