Matlab ODE45 function debugging

%Math462 HW3 Q5b 
fake_data_vals = [86; 86; 117; 120; 130; 135;  169;  179;  224;  230;  242;  234;  244;  245; 
     270;  271;  309; 354;  438;  476;  508;  528;  599;  759;  779;  844;  888;  964; 
     1048;  1093;  1201;  1322;  1437;  1617;  1766;  1835;  1963;  2115;  2225;  2458; 
     2599;  3052;  3944;  4269;  4366;  4963;  5325;  5843;  6242;  6553;  7157;  7470;  8011; 
     8376;  8973;  9191;  9911; 10114; 13676;  13540;  13015;  13241;  14068; 14383;  15113; 
     15319;  15901;  16899;  17111;  17908;  18569;  19463;  20171;  20712;  21261;  21689;  22057;
     22460;  22859;  23218;  23694;  23934;  24247;  24666;  24872;  25178;  25515;  25791;  26044; 
     26277;  26593;  26724;  26933;  27013;  27145;  27237;  27305;  27443;  27514;  27573;  27642;
     27705;  27748;  27862;  27929;  27952;  28005;  28073;  28147;  28220;  28295;  28388;  28421; 
     28454;  28476;  28539;  28571;  28599;  28598;  28601;  28601;  28601;  28604;  28601;  28601; 
     28601;  28601];


fake_data_times  = [1:1:127]';
fake_data = [fake_data_times, fake_data_vals];


%Now we need to run a simulation to compare to the data.
%Initial guesses for the parameters.
%x is our t in the homework problem
alpha0 =0.05;
beta0=0.2;
gamma0=0.125;
N0=50000;
I0=86/N0;
E0=2*I0;
S0=1-I0-E0;
R0=0;
y0=N0*alpha0*E0;
x0=0.2;
%x0 = 0.5;
p0 = [alpha0;beta0;gamma0;I0;E0;S0;R0;y0];
m=[alpha0;beta0;gamma0;N0;I0;E0;S0;R0;y0];
x=[S0,E0,I0,R0,y0];
%do an inital run of the model with the initial values to see how it
%compares to the data
tspan=[1:1:127];
options = odeset('AbsTol',1e-8,'RelTol',1e-8);
df = @(t,x) model_ode(t,x,alpha0,beta0,gamma0);
[~, xsol] = ode45(df, tspan, x0,options);

%plot
clf();
subplot(1,2,1);
plot(tspan,xsol);
hold on;
plot(data_times,data_vals,'o');
title('Initial guess of the model');
xlabel('Time');
ylabel('x(t)');

%now run fminsearch -- maximize the likelihood
minimize_me = @(p) cost_fcn(p,fake_data);
%set some options for fminsearch
options = optimset('Display','iter','MaxFunEvals',5000,'MaxIter',5000);
%run fminsearch
[parvals, NLL] = fminsearch(minimize_me, p0, options);
%spit out the values that best match the data
alphanew = parvals(1);
betanew=parvals(2);
gammanew=parvals(3);
I0new=parvals(4);
E0new=parvals(5);
S0new=parvals(6);
R0new=parvals(7);
y0new=parvals(8);


%display the output of the minimization
disp(parvals);

%re-run the system with the new parameter values
options = odeset('AbsTol',1e-8,'RelTol',1e-8);
df = @(t,x) model_ode(t,x,alphanew,betanew,gammanew);
[~, xsol_new] = ode45(df, tspan, x0_new,options);

%plot the new fit and compare to data
subplot(1,2,2);
plot(tspan,xsol_new);
hold on;
plot(data_times,data_vals,'o');
title('Least Squares Fit');
xlabel('Time');
ylabel('x(t)');

%Below we define two functions that we'll need to do this estimation.
%First, define the differential equation we want to solve, letting it take
%a vector of parameters as an input argument.

function system = model_ode(~,x,params)    
    alpha = params(1); 
    beta=params(2);
    gamma=params(3);
    N=50000;
     
    Sdot=-beta*x(1)*x(2);
    Edot=beta*x(1)*x(2)-alpha*x(3);
    Idot=alpha*x(3)-gamma*x(2);
    Rdot=gamma*x(2);
    ydot=N*alpha*x(3);
    system=[Sdot;Idot;Rdot;Edot;ydot];
end


function NLL = cost_fcn(params, data)    
%define some values    
    alpha = params(1);    
    beta=params(2);
    gamma=params(3);
    N=50000;
    I=params(4);
    E=params(5);
    S=params(6);
    R=params(7);
    y=params(8);
    %interpret the data   
    times = data(:,1);    
    vals = data(:,2);
    m1=[alpha,beta,gamma];
    %find the model values    
    options = odeset('AbsTol',1e-8,'RelTol',1e-8);    
    df = @(t,x) model_ode(t,x,alpha,beta,gamma);    
    [~, xsol] = ode45(df, times, x0,options);
    
    %compare the model to the data    
    %if we aren't using the sum of squares as the log likelihood, we would    
    %change the two lower lines to represent whatever likelihood function    
    %we want to use    
    res = (xsol - vals).^2; 
    %compute the square of the differences    
    NLL = sum(res);         
    %compute the sum of the squares
end

Could someone help me figure out where are the bugs? I think matlab just does not run. I am not sure about how to use ODE45. I might have defined something wrong but I cannot figure it out.

Answer 1

[Cris LaPierre]

There are a couple issues with how you've set up your odefunc. You may find this example from the ode45 documentation page helpful.

1. Your initial conditions need to be the same length as your output (5). Your x0 is a single value.
2. Your function declaration has to match your function call. Your odefunc expects params to be a vector, but you pass in the values as separate inputs (separated by commas). Perhaps you meant to use varargin? It's not necessary, though.

Here's how I would modify just the ode45 relevant code

alpha0 =0.05;
beta0=0.2;
gamma0=0.125;

x=[0.9948, 0.0034, 0.0017, 0, 8.6000];
tspan=[1 127];
options = odeset('AbsTol',1e-8,'RelTol',1e-8);
df = @(t,x) model_ode(t,x,alpha0,beta0,gamma0);
[t, xsol] = ode45(df, tspan, x,options);

plot(t,xsol)
legend

function system = model_ode(~,x,alpha,beta,gamma)    
    N=50000;
     
    Sdot=-beta*x(1)*x(2);
    Edot=beta*x(1)*x(2)-alpha*x(3);
    Idot=alpha*x(3)-gamma*x(2);
    Rdot=gamma*x(2);
    ydot=N*alpha*x(3);
    system=[Sdot;Idot;Rdot;Edot;ydot];
end

[Zhukun Wang]

%Math462 HW3 Q5b

fake_data_vals = [86; 86; 117; 120; 130; 135; 169; 179; 224; 230; 242; 234; 244; 245;

270; 271; 309; 354; 438; 476; 508; 528; 599; 759; 779; 844; 888; 964;

1048; 1093; 1201; 1322; 1437; 1617; 1766; 1835; 1963; 2115; 2225; 2458;

2599; 3052; 3944; 4269; 4366; 4963; 5325; 5843; 6242; 6553; 7157; 7470; 8011;

8376; 8973; 9191; 9911; 10114; 13676; 13540; 13015; 13241; 14068; 14383; 15113;

15319; 15901; 16899; 17111; 17908; 18569; 19463; 20171; 20712; 21261; 21689; 22057;

22460; 22859; 23218; 23694; 23934; 24247; 24666; 24872; 25178; 25515; 25791; 26044;

26277; 26593; 26724; 26933; 27013; 27145; 27237; 27305; 27443; 27514; 27573; 27642;

27705; 27748; 27862; 27929; 27952; 28005; 28073; 28147; 28220; 28295; 28388; 28421;

28454; 28476; 28539; 28571; 28599; 28598; 28601; 28601; 28601; 28604; 28601; 28601;

28601; 28601];

fake_data_times = [1:1:127]';

fake_data = [fake_data_times, fake_data_vals];

%Now we need to run a simulation to compare to the data.

%Initial guesses for the parameters.

alpha0 =0.05;

beta0=0.2;

gamma0=0.125;

N0=50000;

I0=86/N0;

E0=2*I0;

S0=1-I0-E0;

R0=0;

y0=N0*alpha0*E0;

p0 = [alpha0;beta0;gamma0;I0;E0;S0;R0;y0];

x0=[S0,E0,I0,R0,y0];

%do an inital run of the model with the initial values to see how it

%compares to the data

tspan=[1:1:127];

options = odeset('AbsTol',1e-8,'RelTol',1e-8);

%df = @(t,x) model_ode(t,x,alpha0,beta0,gamma0);

%[~, xsol] = ode45(df, tspan, x0,options);

df = @(t,x) model_ode(t,x,alpha0,beta0,gamma0);

[t, xsol] = ode45(df, tspan, x0,options);

%plot

clf();

subplot(1,2,1);

plot(tspan,xsol);

hold on;

plot(data_times,data_vals,'o');

title('Initial guess of the model');

xlabel('Time');

ylabel('x(t)');

%now run fminsearch -- maximize the likelihood

minimize_me = @(p) cost_fcn(p,fake_data);

%set some options for fminsearch

options = optimset('Display','iter','MaxFunEvals',5000,'MaxIter',5000);

%run fminsearch

[parvals, NLL] = fminsearch(minimize_me, p0, options);

%spit out the values that best match the data

alphanew = parvals(1);

betanew=parvals(2);

gammanew=parvals(3);

I0new=parvals(4);

E0new=parvals(5);

S0new=parvals(6);

R0new=parvals(7);

y0new=parvals(8);

x_new=[I0new E0new S0new R0new y0new];

%display the output of the minimization

disp(parvals);

%re-run the system with the new parameter values

options = odeset('AbsTol',1e-8,'RelTol',1e-8);

%df = @(t,x) model_ode(t,x,alphanew,betanew,gammanew);

%[~, xsol_new] = ode45(df, tspan, x0_new,options);

df = @(t,x) model_ode(t,x,alphanew,betanew,gammanew);

[~, xsol_new] = ode45(df, tspan, x0_new,options);

%plot the new fit and compare to data

subplot(1,2,2);

plot(tspan,xsol_new);

hold on;

plot(data_times,data_vals,'o');

title('Least Squares Fit');

xlabel('Time');

ylabel('x(t)');

%Below we define two functions that we'll need to do this estimation.

%First, define the differential equation we want to solve, letting it take

%a vector of parameters as an input argument.

function system = model_ode(~,x,alpha,beta,gamma)

N=50000;

Sdot=-beta*x(1)*x(2);

Edot=beta*x(1)*x(2)-alpha*x(3);

Idot=alpha*x(3)-gamma*x(2);

Rdot=gamma*x(2);

ydot=N*alpha*x(3);

system=[Sdot;Idot;Rdot;Edot;ydot];

end

function NLL = cost_fcn(params, data)

%define some values

alpha = params(1);

beta=params(2);

gamma=params(3);

N=50000;

I=params(4);

E=params(5);

S=params(6);

R=params(7);

y=params(8);

%interpret the data

times = data(:,1);

vals = data(:,2);

%find the model values

options = odeset('AbsTol',1e-8,'RelTol',1e-8);

df = @(t,x) model_ode(t,x,alpha,beta,gamma);

[~, xsol] = ode45(df, times, x0,options);

%compare the model to the data

%if we aren't using the sum of squares as the log likelihood, we would

%change the two lower lines to represent whatever likelihood function

%we want to use

res = (xsol - vals).^2;

%compute the square of the differences

NLL = sum(res);

%compute the sum of the squares

end

THIS IS MY LATEST CODE. I think the code does graph something but it does not display the parameters after I said:

%display the output of the minimization

disp(parvals);

There is no parvals coming out. Just wondering if you could help test where are the bugs? Thanks a lot.

[Cris LaPierre]

Are you getting an error message? If so, please share the complete error message (all the red text).

[Zhukun Wang]

Oh, sorry about that. All data are following: (Note the 1st column is time and 2ed is the data values)

data = [0 86; 1 86; 2 117; 6 120; 7 130; 8 135; 13 169; 16 179; 23 224; 27 230; 29 242; 36 234; 41 244; 50 245; 59 270; 63 271;

64 309; 69 354; 72 438; 77 476; 78 508; 85 528; 91 599; 99 759; 104 779; 105 844; 111 888; 113 964;

118 1048; 121 1093; 125 1201; 128 1322; 131 1437; 132 1617; 136 1766; 140 1835; 141 1963; 143 2115;

147 2225; 149 2458; 150 2599; 156 3052; 165 3944; 167 4269; 171 4366; 175 4963; 177 5325; 181 5843;

183 6242; 185 6553; 190 7157; 192 7470; 197 8011; 199 8376; 204 8973; 206 9191; 211 9911;214 10114;

218 13676; 220 13540; 225 13015; 227 13241; 232 14068; 234 14383; 239 15113; 241 15319; 246 15901;

248 16899; 253 17111; 260 17908; 267 18569; 274 19463; 281 20171; 288 20712; 295 21261; 302 21689;

309 22057; 316 22460; 323 22859; 330 23218; 337 23694; 344 23934; 351 24247; 358 24666; 365 24872;

372 25178; 379 25515; 386 25791; 393 26044; 400 26277; 407 26593; 414 26724; 421 26933; 428 27013;

435 27145; 442 27237; 449 27305; 456 27443; 463 27514; 470 27573; 477 27642; 484 27705; 491 27748;

498 27862; 505 27929; 512 27952; 519 28005; 527 28073; 534 28147; 541 28220; 548 28295; 555 28388;

562 28421; 569 28454; 576 28476; 583 28539; 590 28571; 596 28599; 603 28598; 610 28601; 617 28601;

624 28601; 631 28604; 638 28601; 645 28601; 652 28601; 659 28601];

[Zhukun Wang]

function or variable 'x0' can't be recognized.

error Math462HW3Q5b>cost_fcn (line 124)

[~, xsol] = ode45(df, times, x0,options);

error Math462HW3Q5b>@(p)cost_fcn(p,fake_data) (line 55)

minimize_me = @(p) cost_fcn(p,fake_data);

error fminsearch (line 200)

fv(:,1) = funfcn(x,varargin{:});

error Math462HW3Q5b (line 59 )

[parvals, NLL] = fminsearch(minimize_me, p0, options);

These are my errors.

[Cris LaPierre]

See this page.