- Your initial conditions need to be the same length as your output (5). Your x0 is a single value.
- 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.
Matlab ODE45 function debugging
5 views (last 30 days)
Show older comments
%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.
0 Comments
Answers (1)
Cris LaPierre
on 20 Mar 2021
There are a couple issues with how you've set up your odefunc. You may find this example from the ode45 documentation page helpful.
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
5 Comments
See Also
Categories
Find more on Numbers and Precision in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!