Lsqcurvefit and ode45 to fit data to system of differential equations

Hi! I am currently running into issues with using lsqcurvefit to produce parameters that fit a specific type of disease model. My goal is to fit a system of differential equations to actual data that I have collected while optimizing specific parameters (b, g, d) in the system of differential equations. The main issue I am running into is that my code compiles and runs, but it continues to output my initial gusses for b, g, d rather than different values. I have run into this issue with a few other methods using fminuc() and a meansquarederror function (to fit my data), but nothing appears to work. Any help would be appreciated, thanks!

My code works by extracting the S, I, R, D data (for the disease model) and setting up initial conditions for ode45() to solve the system of differential equations. From here, initial gusses of the the b, g, d parameters are made and lsqcurvefit is called to a function which contains the ode45 and the system of differential equations.

File #1:

% Extracting Data: 
tspan = [1:30]'; % Days 1 to 30
SSS = italyData(1:30,6);
III = italyData(1:30,7);
RRR = italyData(1:30,9);
DDD = italyData(1:30,8);

% Initial Conidtions [S(0) I(0) R(0) D(0)]:
IC = [SSS(1) III(1) RRR(1) DDD(1)];

% Total Data Set
xdata = tspan;
ydata = [SSS(1:30,1) III(1:30,1) RRR(1:30,1) DDD(1:30,1)];

% Initial Guesses (of parameters):
IG = [0.2 0.031 0.048];

%% LSQCURVEFIT
fun = @(x,t) Testv3Solver(tspan, IC, IG);
output = lsqcurvefit(fun, IG, xdata, ydata);
b = output(1)
g = output(2)
d = output(3)

File #2:

function state_out = Testv3Solver(tspan,IC, IG)

fun2 = @(t,IC) odeFun(tspan,IG,IC);
[tout,state_out] = ode45(fun2,tspan, IC);
    
    function state_dot = odeFun(tspan, IG, IC)
    b = IG(1);
    g = IG(2);
    d = IG(3);
    
    state_dot = zeros(4,1);
    S = IC(1);
    I = IC(2);
    R = IC(3);
    D = IC(4);
    N = 60e6;
    
    state_dot(1) = (-b*S*I)/N;
    state_dot(2) = ((b*S*I)/N) - g*I - d*I;
    state_dot(3) = g*I;
    state_dot(4) = d*I;
    end

end

Answer 1

[Star Strider]

Note that neither ‘x’ or ‘t’ are explicity arguments to ‘Testv3Solver’, so nothing is passed to it and nothing ever changes —

fun = @(x,t) Testv3Solver(tspan, IC, IG);

I have posted a series of Answers on this topic. See Parameter Estimation for a System of Differential Equations for a working example.

.

[Wonderboy130]

Thanks for your quick response! I followed the Igor_Moura.m file you provided in response to the question that was posted in the link you provided, but I am still running into trouble with indexing, specifically the ydata in the system of differential equations portion, and I cannot seem to figure out why as the indexing of ydata makes sense ydata is a 30x4 matrix and calling ydata(1,1), ydata(1,2), ... all work).

Error:

File #1:

% Extracting Data: 
tspan = [1:30]'; % Days 1 to 178
SSS = italyData(1:30,6);
III = italyData(1:30,7);
RRR = italyData(1:30,9);
DDD = italyData(1:30,8);

% Initial Conidtions [S(0) I(0) R(0) D(0)]:
IC = [SSS(1) III(1) RRR(1) DDD(1)];

% Total Data Set
xdata = tspan;
ydata = [SSS III RRR DDD];

% Initial Guesses (of parameters):
IG = [0.2 0.031 0.048];

%LSQCURVEFIT
output = lsqcurvefit('Testv3Solver', IG, xdata, ydata);
b = output(1)
g = output(2)
d = output(3)

File #2:

function state_out = Testv3Solver(guess,tspan)

IC = [59999997 3 0 0]; % Initial conditions
[tout,output] = ode45(@odeFun,tspan, IC);
    
    function outHere = odeFun(tspan, ydata) 
    b = guess(1);
    g = guess(2);
    d = guess(3);
    
    state_dot = zeros(4,1);
    S = ydata(1,1); % <--
    I = ydata(1,2); % <--
    R = ydata(1,3); % <--
    D = ydata(1,4); % <--
    N = 60e6;
    
    state_dot(1) = (-b*S*I)/N;
    state_dot(2) = ((b*S*I)/N) - g*I - d*I;
    state_dot(3) = g*I;
    state_dot(4) = d*I;
    outHere = state_dot;
    end
state_out = output;
end

[Star Strider]

Each one of these —

    S = ydata(1,1); % <--
    I = ydata(1,2); % <--
    R = ydata(1,3); % <--
    D = ydata(1,4); % <--

are scalars, not matrices.

Try this instead —

    S = ydata(1); % <--
    I = ydata(2); % <--
    R = ydata(3); % <--
    D = ydata(4); % <--

I cannot run your code to test this, however that should work.

Also, just to be safe, use function handle representation here instead of the single quotes —

    output = lsqcurvefit(@Testv3Solver, IG, xdata, ydata);

.

[Wonderboy130]

Thank you! I think changing that allowed everything to work as I would expect.

[Star Strider]

As always, my pleasure!