How to handle discrete, non-periodic right-hand side of ODE in Matlab

I am trying to migrate an existing Simulink simulation to pure Matlab, which is necessary for several reasons. It is a stiff system of 10 - 50 ordinary, first-order, non-linear differential equations and a couple of algebraic equations, formulated as a level-2 C-Code S-function.

A simplified version of the corresponding Simulink flowsheet (strangely enough, I cannot upload slx files here):

Which produces the following plot:

Consider this MWE for illustration purposes (Matlab version of the Simulink sheet):

% simpleStorageTankModel_Matlab Test of Matlab Solver for Simulink problem
contOutlet = @(t) 0.03*t.^2 - 0.05*t + 0.01;
t1 = 0:1e-2:12;
stairs(t1, periodicDiscreteInlet(t1), 'b'); hold on

% Formulate ODE
solverOptions = odeset('MaxStep', 1e-2);
tic; fillingLevel = ode15s(@(t,y) periodicDiscreteInlet(t) - contOutlet(t),...
    t1, 0.5, solverOptions); toc
plot(t1, 0.03*t1.^2 - 0.05*t1 + 0.01, 'k');
plot(fillingLevel.x, fillingLevel.y, 'r');
legend({'Inlet', 'Outlet', 'Filling level'}, 'Location', 'NorthWest');

function perDisInp = periodicDiscreteInlet(t)
    amplitude = 4;
    phaseDelay = 0;
    pulseWidth = 0.5;
    t_cyc = t - floor(t);
    if numel(t) > 1
        numberOfTimePoints = length(t_cyc);
        perDisInp = zeros(numberOfTimePoints, 1);
        for timeIndex = 1:numberOfTimePoints
            if (t_cyc(timeIndex) >= phaseDelay) && (t_cyc(timeIndex) <...
                1 - pulseWidth + phaseDelay)
                perDisInp(timeIndex) = amplitude;
            end
        end
    else
        if (t_cyc >= phaseDelay) && (t_cyc < 1 - pulseWidth + phaseDelay)
            perDisInp = amplitude;
        else
            perDisInp = 0;
        end
    end
end

An annoying problem is that, if do not set the MaxStep option to something below, say, 1e-2, then the solution beyond t > 6 becomes inaccurate, as if the solver is missing some of the square waves. This seems to have occured to others before .

However, my main problem is with ode15s ignoring the time vector argument I am providing, probably due to the right-hand side of the ODE ( periodicDiscreteInlet) being discrete, i.e., not continuously differentiable (in my simulation, the input is actually a non-periodic square wave). Thus, although t1 is a 1x1201 double vector, I get a 1x1224 double solution. I know I could write a function to extract the values I need, but this would add a significant performance penalty. In my case, Simulink takes just about a tenth of the simulation time of pure Matlab!

Is there a built-in function to "round off" the square-wave edges, making the input continuous (I don't mean as round as via Fourier transform)? Or a stiff Matlab solver that can handle discrete right-hand sides? Why can Simulink handle this problem with such ease, but Matlab can't? Any help is appreciated.

Answer 1

[winkmal]

Since this was inspired by the valuable comment of Walter Roberson to my original answer, I am posting this solution as a separate answer.

So it turns out the for-loop solution actually performs 6-8x faster than interp1 inside the ODE RHS function! Here's a plot of both solutions:

Code below.

% simpleStorageTankModel_Matlab Test of Matlab Solver for Simulink problem
clc; clearvars

p1 = [0.03 -0.05 0.01]; % contOutlet = @(t) 0.03*t.^2 - 0.05*t + 0.01;
t1 = 0:0.5:12;
pdi = periodicDiscreteInlet(t1);
discreteInput = [t1', pdi, 0.5*ones(25,1)];

%% Solve ODE in a loop without interp1
fillingLevel1 = zeros(length(t1), 1);
initialFillingLevel = 0.5;
fillingLevel1(1) = initialFillingLevel;
tic
for timeInd = 1:length(t1)-1
    odeRhs1 = @(t, y) periodicDiscreteInlet(t1(timeInd)) - polyval(p1, t);
    [~, solVec] = ode15s(odeRhs1, linspace(t1(timeInd), t1(timeInd+1), 3), ...
        initialFillingLevel);
    fillingLevel1(timeInd+1) = solVec(end); % Returns 3 values
    initialFillingLevel = solVec(end);
end
toc

%% Solve ODE with interp1
solverOptions = odeset('MaxStep', 1e-2);
modelOdeFun = @(t, y) odeRhs2(t, y, discreteInput, p1);
tic; [~, fillingLevel2] = ode15s(modelOdeFun, t1, 0.5, solverOptions); toc

%% Plot results: V1 as loop, V2 with interpolation
set(groot, 'DefaultLineLineWidth', 1.25);
set(groot, 'DefaultStairLineWidth', 1.25);

figure();
stairs(discreteInput(:,1), discreteInput(:,2), 'b');
hold on; grid on
plot(t1, polyval(p1, t1), 'k');
plot(t1, fillingLevel1, 'r--');
plot(t1, fillingLevel2);
legend({'Inlet', 'Outlet', 'Filling level-loop', 'Filling level-interp1'}, ...
    'Location', 'NorthWest');
title('RHS with discrete jumps, loop vs. interpolation');
set(gca, 'FontSize', 12);

%% Separate functions
function dy = odeRhs2(t, y, discreteInput, p1)
    % Interpolation as specified in Matlab documentation
    % https://mathworks.com/help/matlab/ref/ode45.html#bu3l43b 
    discreteInput_t = interp1(discreteInput(:,1), discreteInput(:, 2), t, 'previous', 0);
    
    dy = discreteInput_t - polyval(p1, t);
end

function perDisInp = periodicDiscreteInlet(t)
    amplitude = 4;
    phaseDelay = 0;
    pulseWidth = 0.5;
    t_cyc = t - floor(t);
    if numel(t) > 1
        numberOfTimePoints = length(t_cyc);
        perDisInp = zeros(numberOfTimePoints, 1);
        for timeIndex = 1:numberOfTimePoints
            if (t_cyc(timeIndex) >= phaseDelay) && (t_cyc(timeIndex) <...
                    1 - pulseWidth + phaseDelay)
                perDisInp(timeIndex) = amplitude;
            end
        end
    else
        if (t_cyc >= phaseDelay) && (t_cyc < 1 - pulseWidth + phaseDelay)
            perDisInp = amplitude;
        else
            perDisInp = 0;
        end
    end
end

Answer 2

[Jan]

Matlab's ODE solvers are designed to handle smooth functions only, see http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. The only reliable and numerically correct solution is to use event functions to stop and restart the integration at all discontinuities.

Answer 3

[winkmal]

After some research, I've found MathWorks official answer, and the corresponding paragraph from ode45 documentation. Note that interp1 does linear interpolation by default, so I changed that to previous.

However, performance in Matlab is much slower than in Simulink.

For completeness, here is my (updated) code, comparing two variants: without and with interpolation, though the former needs an expression f(t) which returns the data at the queried point. Though results are the same (isequal(fillingLevel1, fillingLevel2) gives 1), it takes about twice as long with interpolation, as it takes without it.

The problem with too many return values seemed to be related to the calling syntax of ode15s. It is now fixed.

Very useful comment, Walter Roberson. So you are suggesting a loop to start and stop integration each time a discontinuous jump occurs? This would mean 25 calls to ode15s. Sure this would be faster than using interp1? Maybe I will try soon.

clc; clearvars

p1 = [0.03 -0.05 0.01]; % contOutlet = @(t) 0.03*t.^2 - 0.05*t + 0.01;
t1 = 0:0.5:12;
pdi = periodicDiscreteInlet(t1);
discreteInput = [t1', pdi, 0.5*ones(25,1)];

%% Formulate ODE
solverOptions = odeset('MaxStep', 1e-2);
odeRhs1 = @(t, y) periodicDiscreteInlet(t) - polyval(p1, t);
modelOdeFun = @(t, y) odeRhs2(t, y, discreteInput, p1);
tic; [~, fillingLevel1] = ode15s(odeRhs1, t1, 0.5, solverOptions); toc
tic; [~, fillingLevel2] = ode15s(modelOdeFun, t1, 0.5, solverOptions); toc

%% Plot results: V1 without interpolation
set(groot, 'DefaultLineLineWidth', 1.25);
set(groot, 'DefaultStairLineWidth', 1.25);

figure();
stairs(discreteInput(:,1), discreteInput(:,2), 'b');
hold on; grid on
plot(t1, polyval(p1, t1), 'k');
plot(t1, fillingLevel1, 'r');
legend({'Inlet', 'Outlet', 'Filling level'}, 'Location', 'NorthWest');
title('RHS with discrete jumps, *no* interpolation');
set(gca, 'FontSize', 12);

%% Plot results: V2 with interpolation
figure();
stairs(discreteInput(:,1), discreteInput(:,2), 'b');
hold on; grid on
plot(t1, polyval(p1, t1), 'k');
plot(t1, fillingLevel2, 'r');
legend({'Inlet', 'Outlet', 'Filling level'}, 'Location', 'NorthWest');
title('RHS with discrete jumps, *previous* interpolation');
set(gca, 'FontSize', 12);

%% Separate functions
function dy = odeRhs2(t, y, discreteInput, p1)
    % Interpolation as specified in Matlab documentation
    % https://mathworks.com/help/matlab/ref/ode45.html#bu3l43b 
    discreteInput_t = interp1(discreteInput(:,1), discreteInput(:, 2), t, 'previous', 0);
    
    dy = discreteInput_t - polyval(p1, t);
end

function perDisInp = periodicDiscreteInlet(t)
    amplitude = 4;
    phaseDelay = 0;
    pulseWidth = 0.5;
    t_cyc = t - floor(t);
    if numel(t) > 1
        numberOfTimePoints = length(t_cyc);
        perDisInp = zeros(numberOfTimePoints, 1);
        for timeIndex = 1:numberOfTimePoints
            if (t_cyc(timeIndex) >= phaseDelay) && (t_cyc(timeIndex) <...
                    1 - pulseWidth + phaseDelay)
                perDisInp(timeIndex) = amplitude;
            end
        end
    else
        if (t_cyc >= phaseDelay) && (t_cyc < 1 - pulseWidth + phaseDelay)
            perDisInp = amplitude;
        else
            perDisInp = 0;
        end
    end
end

[Walter Roberson]

Using previous guarantees that the function is not differentiable, so MATLAB will either complain that it has hit a discontinuity or else will have an integration failure at each discontinuity that MATLAB has to back the integration step off a lot until it finds a small enough step size that the error tolerances are satisfied (which can be impossible when yu use previous.

Using linear interpolation or spline can lead one to suspect that there are important discontinuities that are not being accounted for properly, but using 'next' or 'previous' makes it certain that your ode will be broken.

If you plot the values being interpolated with their breakpoints and associated values, and you can visually see any corners or angles in the output, then using that interpolated data will violate the mathematical conditions used by ode45.

... And Yes, we are saying that the support article is wrong.

Any time you have a time based change in value, you should be breaking up the integration into seperate time ranges.