Adding a piecewise time dependent term in system of differential equation

20 views (last 30 days)
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.

Accepted Answer

Thiago Henrique Gomes Lobato
Edited: Thiago Henrique Gomes Lobato on 27 Oct 2019
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
  16 Comments
Walter Roberson
Walter Roberson on 13 Aug 2021
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)'});

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 27 Oct 2019
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.

Categories

Find more on General Applications in Help Center and File Exchange

Products


Release

R2016a

Community Treasure Hunt

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

Start Hunting!