Periodic input to ODE solver

조회 수: 20(최근 30일)
Eric Scott
Eric Scott 2022년 8월 2일
댓글: Eric Scott 2022년 8월 2일
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);
function zr = dInput(t,z)
dCycle = mod(t,1);
if dCycle == 0
input = 1;
input = 0;
zr = input-0.1*z;
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?

채택된 답변

William Rose
William Rose 2022년 8월 2일
편집: William Rose 2022년 8월 2일
Create a smooth pulsatile input signal with a specified non-zero duration. Add it to your ODE. For example, consider the following:
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
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.
  댓글 수: 3
Eric Scott
Eric Scott 2022년 8월 2일
Thank you - you have pointed me in just the direction I needed.

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Walter Roberson
Walter Roberson 2022년 8월 2일
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.
  댓글 수: 1
Eric Scott
Eric Scott 2022년 8월 2일
Understood, thank you for the clarification.

댓글을 달려면 로그인하십시오.





Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by