How to Incorporate a Piecewise Function into an ODE System?
5 views (last 30 days)
Show older comments
The function p(t) in my system of ordinary differential equations is a piecewise function. I have tried two different approaches for incorporating the piecewise function in Matlab, however, none have worked.
For my first approach, I defined p(t) in the ODE file (a nested function):
function dy = ODE_v1(t,y)
dy(1,1) = j(t)*y(1)-k(t)*y(1)-0.25*p(t)*y(1);
dy(2,1) = j(t)*y(1)-k(t)*y(2)-0.50*p(t)*y(2);
dy(3,1) = j(t)*y(2)-k(t)*y(3)-0.75*p(t)*y(3);
% nested functions
function y = j(t)
y = 0.31*sin(t*2*pi/365);
end
function y = k(t)
y = 0.91*sin(t*4*pi/365);
end
function y = p(t)
y = ...
(0*t ) .* (t >= 0 & t < 500) + ...
(0*t+0.50 ) .* (t == 500) + ...
(0*t ) .* (t > 500 & t < 750) + ...
(0*t+0.50 ) .* (t == 750) + ...
(0*t ) .* (t > 750);
end
end
For my second approach, I defined p(t) in the ODE solver.
tf = 1000;
dt = 0.01;
t = 0:dt:tf;
y0 = [300;200;100];
p = @(t) q(t); % p(t) is defined here, where q(t) is defined in a separate file q.m
callODE = @(t,y) ODE_v1(t,y,p); % In the ODE file: function dy = ODE_v1(t,y,p)
[t,y] = ode45(callODE,t,y0);
The file q.m:
function y = q(t)
y = ...
(0*t ) .* (t >= 0 & t < 500) + ...
(0*t+0.50 ) .* (t == 500) + ...
(0*t ) .* (t > 500 & t < 750) + ...
(0*t+0.50 ) .* (t == 750) + ...
(0*t ) .* (t > 750);
end
For both approaches, Matlab considered p(t) as p(t) = 0. Are there other approaches for incorporating a piecewise function? Thank you!
1 Comment
Answers (1)
NING LIU
on 7 Jun 2021
Is your problem solved? I have the same kind of problem...
3 Comments
NING LIU
on 7 Jun 2021
Thanks for your comments. I tried to create the piecewise function with Heaviside, and then integrate it into ode(). Fortunately it works!
Walter Roberson
on 7 Jun 2021
Heaviside() is primarily for symbolic work; you would use dsolve() for that, not one of the ode() functions. However, dsolve() is not very likely to find closed form solutions for piecewise functions.
If you defined a function with heaviside or with using conditional computation along the line of sin(x).*(t>2) and you used one of the ode*() functions, then whatever result you got out was very likely wrong. The functions such as ode45() do not always notice that the computation is meaningless, but the mathematical operations they carry out are only valid under the conditions I described above. You might not have received an error message, but you are unlikely to have received the correct answer.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!