Custom periodic function without counter

Question

I am using ode45 to solve a simple ODE:

function dCdt=u_vent(t,C)
if t>   600 &&  t<= 720
    Q=Q2;               
elseif t>   1320    &&  t<= 1440
    Q=Q2;           
elseif t>   2040    &&  t<= 2160
    Q=Q2;           
elseif t>   2760    &&  t<= 2880
    Q=Q2;               
elseif t>   3480    &&  t<= 3600
    Q=Q2;
else
    Q=Q1;
end

V=100;
C_i=400;
S=100;

dCdt=Q/V*C_i+S/V-Q/V*C(1);
return

I use then to solve:

[t,C]=ode45(@u_vent, [0 1*3600], 400);

I would like to create a periodic function such as the one in the picture for Q, 0<t<3600, without using those if statements... any thoughts?

Answer 1

This is a common type of question. Users often wish to discontinuously change a variable within the integration function used by Matlab's ODE solvers. Unfortunately, this is usually a bad idea. These solvers assume that the differential equations are smooth and continuous. At best your code will be slower. At worst you'll have larger errors or even completely wrong results. Using a stiff solver such as ode15s might help a bit, but it too assumes continuity. As your system as specified has analytic solution, the easiest way so "simulate" it might actually be to pass a pulse train trough this function of time.

The case when a parameter changes discontinuously in time, i.e., with respect to the independent variable, is easy to solve. Simply integrate over each time span for the number of periods you want:

t1 = 600;
t2 = 120;
TSpan = [0 t1]; % Initial time vector
NPeriods = 5;   % Number of periods
C0 = 400;       % Initial condition
Q1 = ...
Q2 = ...
V = 100;
C_i = 400;
S = 100;
u_vent = @(t,C,Q)(Q*(C_i-C(1))+S)/V; % Integration function as anonymous function
Cout = C0;       % Array to store solution states
tout = Tspan(1); % Array to store solution times
for i = 1:NPeriods
    [t,C] = ode45(@(t,C)u_vent(t,C,Q1), TSpan, C0);
    tout = [tout;t(2:end)];   % Append time, 2:end used to avoid avoid duplicates
    Cout = [Cout;C(2:end,:)]; % Append state
    TSpan = [0 t2]+t(end);    % New time vector
    C0 = C(end);              % New initial conditions
    [t,C] = ode45(@(t,C)u_vent(t,C,Q2), TSpan, C0);
    tout = [tout;t(2:end)];
    Cout = [Cout;C(2:end,:)];
    TSpan = [0 t1]+t(end);
    C0 = C(end);
end

This allows ode45 to integrate an ODE for a specified time from a set of initial conditions. Then the integration is restarted at the end time of the previous integration with new discontinuous initial conditions or different parameters. This results in a transient. You may not like the look of the code, but this is how it's done. It's even trickier if parameters change discontinuously with respect to state variables.

Optional. Each time you call/restart ode45 the function must figure out an initial step size to use. This is the primary (minimal) cost/overhead of restarting the solver. In some cases you may be able to help out the solver by specifying an 'InitialStep' option based on the last step size used from the previous call. See the ballode demo for further details on this by typing edit ballode in your command window.

A note. The tout and Cout arrays are appended to themselves after each integration step. This is effectively dynamic memory allocation. As long as NPeriods is reasonably small, this likely won't be a problem as dynamic memory allocation in recent versions of Matlab can be very fast and you're only reallocating a few time in large blocks. If you specify fixed step size output (e.g., TSpan = 0:1e-2:600;) then you'll know exactly how much memory to preallocate to tout and Cout.

Answer 2

If you have the Signal Processing Toolbox, you can use something like:

>> t = 0:3600;
>> y = pulstran(t,[660:720:3600],'rectpuls',120);
>> plot(t,y)
>> ylim([-0.1 1.1])

which gives the following (in Octave, should be the same in MATLAB):

You then need to scale y to be between Q1 and Q2 instead of 0 and 1.

Answer 3

Not necessarily the best approach, because continuity assumptions remain broken, but the way to generate a rectangular pulse train without the if chain is:

Q = Q2 + (Q1 - Q2) * (mod(t, period) < t_rise);

which works in both scalar and vector context.