Azzera filtri
Azzera filtri

Info

Questa domanda è chiusa. Riaprila per modificarla o per rispondere.

How to simulate an ODE when one of the parameters is a function of time

1 visualizzazione (ultimi 30 giorni)
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);

Risposte (1)

darova
darova il 4 Giu 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 Commenti
Rose
Rose il 4 Giu 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 il 5 Giu 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

Questa domanda è chiusa.

Tag

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by