Info
This question is closed. Reopen it to edit or answer.
How to simulate an ODE when one of the parameters is a function of time
1 view (last 30 days)
Show older comments
I am simulation a system of ODEs in which one parameter is a function of time (similar to a switch function) and others are constants. Below is the MWE and it gives me various errors.
function dxdt = Eq(~,x,par,h)
t2 = 10;
counter =1;
for t=0:1:20
if t>0 && t<t2
h(:,counter) = 1;
else
h(:,counter) = 0;
end
counter = counter+1;
end
disp(size(h))
dxdt = NaN(3,1);
S = x(1);
I = x(2);
V = x(3);
% creating a vector of parameters
K1 = par(1);
k2 = par(2);
k3 = par(3);
k4 = par(4);
U = par(5);
dxdt(1) = -(h*U*K1+k2)*S; %S
dxdt(2) = h*U*K1*S - k3*I; %I
dxdt(3) = k3*I - k4*V;% V
end
Here is the solver file:
clear all
close all
Tfinal = 20;
dt = 1;
t = 0:dt:Tfinal;
disp(size(t))
K1 = 0.5;
k2 = 0.1;
k3 = 0.1;
k4 = 0.05;
U = 1;
t1 = 1;
t2 = 10;
counter =1;
for t=0:1:20
if t>0 && t<t2
h(:,counter) = 1;
else
h(:,counter) = 0;
end
counter = counter+1;
end
disp(size(h))
par = [k1 k2 k3 k4 U];
S0 = 1;
I0 = 0;
V0 = 0;
x0 = [S0 I0 V0]; % you also have to add it to the IC vector
sol = ode23(@Eq,[0 Tfinal],x0,par,h); % solving the model
% size(sol)
x1 = deval(sol,t);
x = x1';
% size(x)
S = x(:,1);
I = x(:,2);
V = x(:,3);
0 Comments
Answers (1)
darova
on 4 Jun 2019
Read how to Pass Extra Parameters in help ode23:
sol = ode23(@(t,x)Eq(t,x,par),[0 Tfinal],x0); % solving the model
Passing trash parameter?
function dxdt = Eq(~,x,par,h) % why do you pass 'h' parameter if you declare it in the function
Why this part of code is repeated in main code? I suggest you to remove it from main code and change it in the function
% for t=0:1:20
% if t>0 && t<t2
% h(:,counter) = 1;
% else
% h(:,counter) = 0;
% end
% counter = counter+1;
% end
% replace with this
h = t>0 && t<t2;
Also read about function handle. If you have simple function you can write it in main code
Eq = @(t,x) [-((0<t&&t<t2)*U*K1+k2)*x(1); ...%S
(0<t&&t<t2)*U*K1*x(1)-k3*x(2); ...%I
k3*x(2) - k4*x(3)];% V
sol = ode23(@(t,x)Eq(t,x,par),[0 Tfinal],x0); % solving the model
2 Comments
Rose
on 4 Jun 2019
I have a similar model and same situation. I almost understoof all your points but how does
h = t>0 && t<t2;
give a step function h(t) which is 1 when 0<t<10 and otherwise 0?
darova
on 5 Jun 2019
't' variable is a scalar in a function. Imagine that Eq function repeats and 't' variable changes every time so 'h' variable would be different
This question is closed.
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!