Having Trouble Using For-Loop Over Distributed Range with ODE45

3 次查看（过去 30 天）
Thomas Veith
Thomas Veith 2019-7-16
评论： Thomas Veith 2019-7-17
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!

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2019-7-17
编辑：Torsten 2019-7-17
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 个评论
显示 3更早的评论
Torsten
Torsten 2019-7-17
编辑：Torsten 2019-7-17
delta1(t) returns true (=1) if t<=1200, false (=0) else,
delta0(t) returns true (=1) if t > 1200, false (=0) else.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by