Clear Filters
Clear Filters

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)
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);

Answers (1)

darova
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
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
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.

Tags

Products

Community Treasure Hunt

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

Start Hunting!