Conditional with time in ode45
14 views (last 30 days)
Show older comments
Brilliant Purnawan
on 10 Jan 2021
Answered: Steven Lord
on 10 Jan 2021
Hi, how would you change a variable in a function when 't' in ODE45 reaches a certain value? I have already tried formatting it like the code below, but that doesn't work I have also tried using the if statement in the function but then it simply won't read the 't' value. So if I want to change upsilon after a certain time has passed in ODE45 how would I approach this?
clear all, close all, clc
mu = 6.25e10^-3;
eta = 6.25e10^-3;
beta = 0.03813350836;
gamma = 0.01165323939;
delta = 1/7;
alpha = 0.00513392178;
tspan = [1:1:500];
y0 = [0.99 0 0.01 0 0];
[t,y] = ode45(@(t,y)odefcn(y,mu,eta,upsilon,beta,gamma,delta,alpha),tspan,y0);
if t>300
upsilon = 0.001;
else
upsilon = 0;
end
plot(t,y)
legend('S','E','I','R','D')
function dydt = odefcn(y,mu,eta,upsilon,beta,gamma,delta,alpha)
dydt = zeros(5,1);
dydt(1) = mu - y(1)*(beta*y(3)+mu+upsilon);
dydt(2) = beta*y(1)*y(3)-y(2)*(mu+delta);
dydt(3) = delta*y(2)-y(3)*(mu+gamma+alpha);
dydt(4) = gamma*y(3)+upsilon*y(1)-eta*y(4);
dydt(5) = alpha*y(3);
end
0 Comments
Accepted Answer
Alan Stevens
on 10 Jan 2021
Include upsilon within the odefcn
mu = 6.25e10^-3;
eta = 6.25e10^-3;
beta = 0.03813350836;
gamma = 0.01165323939;
delta = 1/7;
alpha = 0.00513392178;
tspan = [1:1:500];
y0 = [0.99 0 0.01 0 0];
[t,y] = ode45(@(t,y)odefcn(t, y,mu,eta,beta,gamma,delta,alpha),tspan,y0);
plot(t,y)
legend('S','E','I','R','D')
function dydt = odefcn(t,y,mu,eta,beta,gamma,delta,alpha)
if t>300
upsilon = 0.001;
else
upsilon = 0;
end
dydt = zeros(5,1);
dydt(1) = mu - y(1)*(beta*y(3)+mu+upsilon);
dydt(2) = beta*y(1)*y(3)-y(2)*(mu+delta);
dydt(3) = delta*y(2)-y(3)*(mu+gamma+alpha);
dydt(4) = gamma*y(3)+upsilon*y(1)-eta*y(4);
dydt(5) = alpha*y(3);
end
More Answers (1)
Steven Lord
on 10 Jan 2021
Solve the system from time t = 0 to time t = 300 with the first value of upsilon. Use the value of the solution at t = 300 as the initial condition for a second call to the ODE solver with the updated value of upsilon.
f = @(t, y, k) k*y;
[t1, y1] = ode45(@(t, y) f(t, y, 2), [0 2], 1);
[t2, y2] = ode45(@(t, y) f(t, y, 1), [2 4], y1(end));
plot(t1, y1, 'r-', t2, y2, 'k--')
If you don't have a fixed time when the value of the parameter should change (you instead have a condition you can check to decide when the parameter should change) use an Events function to detect the change. See the ballode example for an illustration of this technique.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!