Periodic input to ODE solver

I would like to run an ODE solver over a time span but provide an additional jump in initial condition on a periodic basis (every day for 112 days). I attempted to do this by using mod(t,1) and setting the time span to be 0:0.1:112, so that every ten steps I could add a static amount to the ODE, and in between it would just have a simple decay term.

t = [0.1:.1:112];

z0 = 0.1;

[t,z] = ode23s(@dInput,t,z0);

plot(t,z);

function zr = dInput(t,z)

dCycle = mod(t,1);
if dCycle == 0
    input = 1;
else
    input = 0;
end

zr = input-0.1*z;

end

mod(t,1) is producing irregular results, sometimes failing to set input for long stretches and other times setting it too frequently. And even when it is set to 1, the output for the equation does not change accordingly.

Is it possible to provide a periodic input to this kind of solver?

Answer 1

[William Rose]

@Eric Scott,

Create a smooth pulsatile input signal with a specified non-zero duration. Add it to your ODE. For example, consider the following:

t=0:.001:6;
tdur=0.1;       %pulse duration (s)
tPer=1.0;       %pulse period, should exceed tdur (s)
tau=0.005;      %pulse time constant
%If tPer=10 and t=0:30, then tmod=[0:5,-5:5,-5:5,-5:0]
tmod=mod(t-tPer/2,tPer)-tPer/2;  %modulo time
%ydown=1./(1+exp((t-tdur)/tau));
yup=1./(1+exp(-tmod/tau));   %this part goes up, from 0 to 1, as tmod crosses 0
ydown=1./(1+exp((tmod-tdur)/tau)); %this part goes down, from 1 to 0, as tmod crosses tdur
input=yup.*ydown;    %this is the pulsatile input

Plot the input versus time:

subplot(2,1,1), plot(t,input,'-r.');
subplot(2,1,2), plot(t,input,'-r.'); xlim([0.95,1.15]); grid on

Plot above shows the signal and a zoomed-in section.

If you want sharper edges, make tau smaller. For a narrower pulse, reduce tdur.

[Eric Scott]

Thank you!

[William Rose]

@Eric Scott,

[edited: correcting typos]

You don't need to use a stiff solver like ode23s. ode45 is fine.

A problem that can occur is that the integrator, when integrating such a simple system - exponential decay at a constant rate - can take large step sizes - so large that it jumps right past the pulse. This can appear to occur almost randomly, producing perplexing results. The cure is to use odeset() to set the options for the integration, specifically to make the maximum step size shorter than the pulse duration. I chose to make the max step size one quarter of the pulse duration:

options=odeset('MaxSize',tDur/4);

To do the above, tDur must be known in the main program. Otherwise, I could have defined tDur in function zr(). Since tDur is in the main program, I need to pass the value of tDur from the main program to function zr(). To do this, I follow the suggestion in ode45() help, under "Pass Extra Parameters...". See my code below and see the example in the ode45 help.

The attached code uses 1 second long pulses (tDur=1.0) that repeat every 5 seconds (tPer=5.0). Time span=[0,40]. The time constant for decay of z (tauZ=10) is 10 seconds, which is consistent with your code (zr=-0.1*z=-z/tauZ). This combination of parameters produces this plot:

which looks a lot like a plot of muscle tension in response to stimulation at a sub-tetanic frequency. This is sometimes called incomplete, or unfused, tetanus - see left hand plot here. (I and many others of a certain age did this experiment with frog muscles, as students. It is no longer common for undergrads to do experiments with actual frog muscles. Simulations are used more now.)

If you make the period (tPer) long compared to tauZ, you get a different looking plot, because z decays almost completely before the next pulse happens. In the example below, tPer=10 and tauZ=2:

Good luck.

[Eric Scott]

Thank you - you have pointed me in just the direction I needed.

Answer 2

[Walter Roberson]

No, that violates the mathematics that the solver uses. You need to stop the function every time there is a discontinuity and restart it. And remember that the solver does not take fixed time steps so do not assume that it will evaluate exactly at the integers.

In your case the discontinuity is predictable just from the time, so loop running first 0.1 to 1, then adjust the boundary conditions then 1.0 to 2.0, adjust conditions, then 2.0 to 3.0... and so on.

Note that when you pass in a vector tspan with more than two elements then the ode*() functions continue to evaluate at whatever times they feel is appropriate to meet the accuracy, but they will only report the results at the time entries you provide.

[Eric Scott]

Understood, thank you for the clarification.