Having Trouble Using For-Loop Over Distributed Range with ODE45

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!

Answer 1

[Torsten]

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.

[Thomas Veith]

Thank you! Now, where do I actually put in my delta values here? I see that you're making it a function of t, which makes sense, but where are the actual values?

[Torsten]

In the two preceeding lines I posted.

delta0 and delta1 must be defined as function handles.

[Thomas Veith]

I'm sorry, I'm confused? All I see is

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

I don't see where you say, for this range, delta1 should equal 1.

[Torsten]

delta1(t) returns true (=1) if t<=1200, false (=0) else,

delta0(t) returns true (=1) if t > 1200, false (=0) else.

[Thomas Veith]

Absolute king