‘What am I doing wrong?’
First, not posting the actual problem, that being:
Warning: Failure at t=7.149052e+00. Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (1.421085e-14) at time t.
so the size discrepancy is due to ode45 not being able to estimate the values for all the time elements.
Second, setting ‘c0’ to have only 3 elements when 5 are required. Fixing two of the parameters within the objective function is not appropriate. A better approach is to constrain them using the parameter bounds vectors.
Note that ‘c0’ needs to have realistic values. All the elements of ‘Sdata’ are between 1 and 17, so the last 3 elements of ‘c0’ in the slightly revised code below needs to reflect that.
I changed the code slightly to reflect this, and also to estimate the initial conditions as part of the other parameters. I leave you to experiment with it:
Time = [0 1 2 3 4 5 6 7 8 9];
Sdata = [3 17 9 4 8 2 1 9 2 4];
% Sdata = Sdata';
figure
plot(Time, Sdata)
grid
c0 = rand(8,1)*10;
[c,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Fit,c0,Time,Sdata);
function S = Fit(c, t, ~)
y0 = c(6:end);
[T,Sv] = ode45(@DifEq, t, y0);
function dS = DifEq(t, y)
ydot = zeros(3,1);
ydot(1) = (c(1)*y(1).^2) + (c(2)*(y(2))) - (c(3)*y(3));
ydot(2) = (c(3)*y(3)) + (c(4)*y(1));
ydot(3) = (c(1)*y(1).^2) - (c(5)*y(2)) + (c(3)*y(3));
dS = ydot;
end
S = Sv(:,3);
end
The problem appears to be that the differential equations may not correctly describe the system, or that the initial parameter estimates are inappropirate. Also, the system appears to be stiff, so ode15s would be an appropriate solver.
Estimating the parameters using the Global Optimization Toolbox ga (genetic algorithm) function would likely provide the best results, since lsqcurvefit requires that the initial parameter estimates are reasonably close to the actual optimal values.
