Adding a piecewise time dependent term in system of differential equation

Problem Statement: I have a system of differential equation and out of which few has time dependent coefficients and in my case i have piecewise time dependence. I am unable to add the time dependence and solve the equation.

Simplified Differential Equation:

My time span runs from 0 to 360 days and my aim is to have certain nonzero value for for half a day and zero for other half.

My unsuccessful attempt:

Function definition: MWE_fn.m

function rk1 = MWE_fn(t,y)
r2=0.5;b2=1;d_A=0.05; c=0.5;
rk1(1)=r2*y(1)*(1-b2*y(1))-c*y(1)*y(2);

rk1(2)= A0(t)- d_A*y(2);

rk1=rk1(:);
end

function fa=A0(t)

    if rem(t,1)==0
        fa=0.04
    else 
        fa=0
    end

end

Function Body: MWE_body.m

timerange= 0:0.5:360;
IC= [0.1,0];%initial conditions 
[t,y] =ode45(@(t,y) MWE_fn(t,y),timerange, IC);

plot(t,y(:,1),'b');
hold on
plot(t,y(:,2),'r');

I was aiming to write some kind of conditional statement for the function to take those values but that didn't go correct. The thing i had in mind was that for half a day we have time at integer values and for other half i have time at 0.5 , so i used a rem function to check that, but i was wrong and it didn't work.

Anyhelp would be appreciated.

Thank You.

Answer 1

[Thiago Henrique Gomes Lobato]

The conditional statement idea is right, the problem was that you considered that the ode45 would evaluate your function only at the time steps you gave, which is not right. To perform the integration many other time intervals are needed, so you just have to adjust your conditional:

function fa=A0(t)
    if t>round(t)
        fa=0.04;
    else 
        fa=0;
    end
    
end

In this case, for every value x.y where y<5, fa = 0.04, which represents basically all your half days

[Vira Roy]

Thanks Thiago. It seems your solution to my original problem is working , so i will accept your solution. I just have a follow up querry.

Now suppose i want to have the nonzero value for n days and zero for other n days and so on. So specifically this would be like nonzero value for first 10 days then zero for next 10 days and so on. How should the conditional statement modify? Will i have to use a for loop also because if statement alone creates problems like this question here where only else statement is considered.

Thank you.

[Thiago Henrique Gomes Lobato]

In your case you can have total information of the actual day based only in your t value, so something like this would work:

function fa=A0(t)
    if t<10
        fa=0.04;
    elseif 10<=t & t<30
        fa=0.08;
    % And so on...
    else
        fa=0;
    end
    
end

You indeed, as Walter Roberson said, may get problems when you have a discontinuity in the same iteration step, this will largely depends about how smooth is your function and how fast do you change the parameter, I had overseen this fact and, ironically, your initial case where there's a change everyday may be strongely more affect for it. The solution of ensuring that for every integration step the condition stays the same can be done by evaluating only partial integrations where the interval will remain in a specific case, as example for your first case of half day on/half day off:

timerange= 0:0.5:360;
IC= [0.1,0];%initial conditions 
for idx=1:length(timerange)-1
    % Integrate
    [tTemp,yTemp] =ode45(@(t,y) MWE_fn(t,y),timerange(idx:idx+1), IC);
    % Save values at the time steps you want to
    t(idx) = tTemp(end); 
    y(idx,:) = yTemp(end,:);
    % Change the actual initial Condition
    IC = y(idx,:);
end
plot(t,y(:,1),'b');
hold on
plot(t,y(:,2),'r');

This gives indeed a smoother results in comparison to the full integration.

[Walter Roberson]

elseif 10<=t<30

would need to be

elseif 10<=t & t<30

[Thiago Henrique Gomes Lobato]

You're right, thanks!

[Gottumukkala Sai Bhavani]

is there a way to plot between this A0(t) and time for the purpose of analysis

[Walter Roberson]

Rewrite A0 to

function fa = A0(t)
    fa = zeros(size(t));
    fa(rem(t,1)==0) = 0.04;
end

Now after your call

[t,y] = ode45(@(t,y) MWE_fn(t,y),timerange, IC);

you can do

a0 = A0(t);
plot(t, a0)

However, you are unlikely to see what you want. Using discontinuous functions like this will give incorrect results. See https://www.mathworks.com/matlabcentral/answers/487643-adding-a-piecewise-time-dependent-term-in-system-of-differential-equation#answer_398437

[Gottumukkala Sai Bhavani]

hey walter, Thank you for response.

May be i should show what i am doing.

for the function shown here i am willing to plot s1(t) vs t to show the inputs on a graph and compare it against my results.

this is part of my main function file FFL_MM

function fa = s1(t)
if t<5
    fa = 0;
elseif t>=5 && t<25
    fa = 1;
elseif t>=60 && t<=80
    fa = 1;
else 
        fa = 0;

end
end

and in my body file FFL_MM_SOL.m

this is being called

T = [0 100];
z0 = [0,0,0];
[t,z] = ode45(@FFL_MM,T,z0);

figure(1);
subplot(3,1,1)
plot(t,z(:,1),'k');
xlim([0 50])
xlabel('Time')

please help me ... i am new to matlab

[Walter Roberson]

function fa = s1(t)
   fa = (t >= 5 & t <25) | (t >= 60 & t <= 80)
end

then after you call ode45, getting back t and z, then

s1vals = s1(t);
plot(t, z(:,1), 'k', t, s1, '--b');

[Gottumukkala Sai Bhavani]

s1vals = s1(t);

i tried this ...

the error showing up is Unrecognized function or variable 's1'.

i might add both the function fa = s1(t); and the ODE call are in twi different files

[Gottumukkala Sai Bhavani]

here are my files

FFL_MM.m

function f =  FFL_MM(t, z)
% omitted parameters 

f(1) = a*s1(t) - kd1*z(1); 
f(2) = b*s2(t) + (k1*z(1)^n1)/(km1^n1 + z(1)^n1) - kd2*z(2); 
f(3) = (k2*z(2)^n2)/(km2^n2 + z(2)^n2) + (k3*z(1)^n3)/(km3^n3 + z(1)^n3) - kd3*z(3);

f = f'; 
end

function fa = s1(t)
if (t >= 5 && t <25) || (t >= 60 && t <= 80)
    fa = 1;
else 
    fa = 0;
end
end

function fb = s2(t)
if (t >= 5 && t <25) || (t >= 60 && t <= 80)
    fb = 1;
else 
    fb = 0;
end
end

and the ode file FFL_MM_sol.m

clear all; close all; clc ;

T = [0 100];
z0 = [0,0,0];
[t,z] = ode45(@FFL_MM,T,z0);


figure(1);
subplot(3,1,1)
plot(t,z(:,1),'k');
xlim([0 50])

and would like to plot (t,s1) and (t,s2) in a seperate subplots and these graphs would look like square/rectangular pulses

[Gottumukkala Sai Bhavani]

seems the following using fplot (changes made im FFL_MM.m) has solved my issue but the following warning message pops up and takes up lot of time for the process to finish

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments

function fa = s1(t)
    fa = zeros(size(t));
if (t >= 5 && t <25) || (t >= 60 && t <= 80)
    fa = 1;
else 
    fa = 0;
end
end
figure(2);
fplot(@(t) s1(t), [0 50]);
ylim([0 2])

[Walter Roberson]

I already posted a version of s1 that would not have that difficulty; https://www.mathworks.com/matlabcentral/answers/487643-adding-a-piecewise-time-dependent-term-in-system-of-differential-equation#comment_1524778

[Muhammad Sinan]

I want to solve a differentia equation with time dependent peacewise defined initial conditions using any ode solver, my main task is to learn how to define any time dependent peacewise defined initial condition that works with ode solve that may or may not be pre-defined in matlab. Like, if we use ode45 or we use our own defined function of Euler scheme.

Hope someone got the idea!

[Walter Roberson]

@Muhammad Sinan

If your piecewise depends only on time and not on state, then call the ode solver repeatedly, each time passing in only a time segment that the piecewise definition would be consistent for. Take the last output of each call, modify it if necessary (e.g., sometimes you want to give an impulse of energy every so-many seconds), and use the result as the boundary conditions for the new call with the next time stretch. In this situation, no event function is needed.

If your piecewise depends upon state, such as "object hit the floor", then you need to use event functions to signal that the condition has been reached. The ode function will terminate at that point. You then take the last output, modify if necessary (e.g., reverse velocity for a bounce), and take the last output time as the new starting time for the next call. Repeat this in a loop until you get to the end of your desired time span.

Do not use piecewise() or if/elseif inside the ode*() functions, and do not use interp1() either (unless you use spline method)... not unless you can guarantee that the first two derivatives of the code you have written would be continuous.

[Muhammad Sinan]

for example, the ode is dxdt = -x; and

if 1<t & t<=1.6

x0 = 0;

elseif 1.6<t & t<=3.2

.

.

end

then how will this condition defined?

[Walter Roberson]

Using a modified form to illustrate holding over values:

x01_1 = 0;
x02_1 = .1;
x01_2 = pi/2;
x01 = [x01_1, x02_1];
[t1, x1] = ode45(@(t,x) [x(2)^2-x(1); x(2)-1/20], [1 1.6], x01);

x02 = x1(end,:); x02(1) = x02_1;
[t2, x2] = ode45(@(t,x) [x(2)^2-x(1); x(2)-1/20], [1.6 3.2], x02);

t = [t1;t2];
x = [x1; x2];
plot(t, x); legend({'x(1)', 'x(2)'});

Answer 2

[Walter Roberson]

You should examine the odeball example to see how to set up events to terminate integration every time there is a change that is not continuous differentiable at least twice. ode45 and related are for continuous systems only.

Or in your case where the determination can be done based just on time, you can work it without events by looping over time intervals.

When you use conditional statements with ode45 and kin then you need to arrange so that the condition is always true or always false within the boundaries of the current integration.

[Vira Roy]

Well, very special thanks to Walter and Thiago for helping me out.