Having Trouble Using For-Loop Over Distributed Range with ODE45

1 view (last 30 days)
Hi all,
So I'm trying to solve an ODE where the values for two of my inputs change over time.
For the first half of my tspan, delta1 should equal 1, and delta0 should equal 0. Then, in the second half, this gets reversed (i.e., delta1=0 and delta0=1). However, the way I've written my code, it only takes the second set of values. That is, it's solving the ODE as is delta1=0 and delta0=0 for the entire time - and I'm not sure why?
Here is my code:
% assign number of parameter sets and create parameters vector
n = 1000;
parameters = zeros(n,6);
% initialize parameters
for k=1:n
alpha0=rand;psa0=29.00;gamma=rand;psi=0.00;beta=rand;mse=0;
tspan=[0 2272];
parameters(k,:) = [alpha0,psa0,gamma,psi*gamma,beta,mse];
end
% solve ODE
for k = 1:size(parameters,1)
options=odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
tspan=[0 2272];
for t=drange(0:1200)
delta1=1;delta0=0;
end
for t=drange(1201:2272)
delta1=0;delta0=1;
end
ODE = @(t,x,delta1,delta0,gamma,psi,beta) [delta1*-gamma*x(1) + delta0*psi*x(1); beta*x(2) - delta1* x(1)*x(2)];
sol = ode45(@(t,x) ODE(t,x,delta1,delta0,parameters(k,3),parameters(k,4),parameters(k,5)), tspan, [parameters(k,1) parameters(k,2)], options);
% calculate MSE and update parameters matrix
a = deval(sol,patientdata1(:,1));
a = a';
cost = (a(:,2)-patientdata1(:,2))./patientdata1(:,2);
mse = cost'*cost;
if min(a(:,2))<= 0.3
mse = mse * 100;
end
parameters(k,6) = mse;
end
Any help would be appreciated. Thanks in advance!

Accepted Answer

Torsten
Torsten on 17 Jul 2019
Edited: Torsten on 17 Jul 2019
delta1 = @(t) (t<=1200)
delta0 = @(t) (t>1200)
ODE = @(t,x,gamma,psi,beta) [delta1(t)*(-gamma)*x(1) + delta0(t)*psi*x(1); beta*x(2) - delta1(t)* x(1)*x(2)];
But a better solution would be to integrate from t=0 to t=1200, store the solution, and call ODE45 for a second time from t=1200 to t=2272 with delta0 and delta1 interchanged.
  5 Comments
Torsten
Torsten on 17 Jul 2019
Edited: Torsten on 17 Jul 2019
delta1(t) returns true (=1) if t<=1200, false (=0) else,
delta0(t) returns true (=1) if t > 1200, false (=0) else.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!