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?

enter image description here

If you have Simulink, then it's easily done in a graphical manner — am304, Dec 06 '13 at 17:22
I assume that you know that your system has an analytic solution? If you're not good with differential equations, you can use Matlab's [`dsolve`](http://www.mathworks.com/help/symbolic/dsolve.html) function: `syms Q V C_i S t C(t) C0;` `dsolve(diff(C,t)==(Q*(C_i-C)+S)/V,C(0)==C0)`. — horchler, Dec 06 '13 at 17:49

horchler · Accepted Answer · 2013-12-06T23:05:50.277

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.

Ah thats fantastic! Thank you very much! It seems like the only sensible solution. If Q were a function, eg the one suggested by @am304 would it still be feasible to so the same thing? — HCAI, Dec 06 '13 at 22:29
If it's a smoothly varying input (e.g., sinusoidal forcing) you may be able to include it in the integration function. Non-smooth inputs like your's may or may not work. By "may work," I mean get results that aren't entirely wrong but still have more error. It depends on the particular ODE and the solver and the solver's settings. Some ODEs need to be carefully integrated when subject to transients or errors will explode. Regardless, errors will be larger than if you use this answer to integrate systems with discontinuities (for the same solver tolerances). And it's a bad habit to get into. — horchler, Dec 06 '13 at 23:01

score 1 · Answer 2 · answered Dec 06 '13 at 17:08

1

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

enter image description here

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

answered Dec 06 '13 at 17:08

am304

13,758
2
22
40

horchler's suggestion is probably better than mine. – am304 Dec 06 '13 at 17:11

score 1 · Answer 3 · answered Dec 06 '13 at 17:32

1

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.

answered Dec 06 '13 at 17:32

Ben Voigt

277,958
43
419
720

Custom periodic function without counter

3 Answers3

Linked