Conditional with time in ode45

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

Accepted Answer

Alan Stevens
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
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.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!