lsqcurvefit (with system of ODEs) does not iterate, it only displays the initial guess value

I have a small problem with LSQCURVEFIT, I think it is easy to solve but for some reason I can't find the solution. I adapted the code below from starstrider "igor-moura.m" file.

The problem is that lsqcurvefit does not start an iteration procedure, but instead only displays the results for the initial guess fitting parameter (theta0). The goal is of course that it changes this parameter to fit the data. In my code, i know the guessing parameter final value "theta" should be close to 41.025 (although i saw 45 gives a better fit). Another old thread suggested to set the options, which I did, but it didn't help.

I use Matlab 2017a version.

function Aspirin_kinetics

function C=kinetics(theta,t)
c0=[0.000996388;	0.002916388]; %sa aa
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(2,1);
    dcdt(1)=-theta(1).*c(1).*c(2);
    dcdt(2)=-theta(1).*c(1).*c(2);
    dC=dcdt;
    end
C=Cv;
end

t = [1
    2
    3
    5
    7
    11
    15
    30];

c = [0.000996388
0.000867473
0.000778297
0.000607295
0.00047941
0.000314246
0.000223062
6.98718E-05]; %SA
c(1:end,2) = 0.002916388 - (0.000996388-c); %AA

theta0=41.025*2;
options = optimset('FinDiffRelStep',1e-3,'DiffMinChange',1e-3,'DiffMaxChange',1);
[theta]=lsqcurvefit(@kinetics,theta0,t,c,[],[],options);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'Location','N')

end

Answer 1

[Star Strider]

It is generally good practice to estimate the initial conditions as well, so I added that (not part of my original code, updated later) and updated the plot so that the colours match (also in my updated versions).

Try this —

Aspirin_kinetics

theta0 = 3×1

82.0500 0.0010 0.0029

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

	Rate Constants:
		Theta(1) = 72.05000
		Theta(2) =  0.00098
		Theta(3) =  0.00299

function Aspirin_kinetics

function C=kinetics(theta,t)
% c0=[0.000996388;	0.002916388]; %sa aa
c0 = theta(2:end);
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(2,1);
    dcdt(1)=-theta(1).*c(1).*c(2);
    dcdt(2)=-theta(1).*c(1).*c(2);
    dC=dcdt;
    end
C=Cv;
end

t = [1
    2
    3
    5
    7
    11
    15
    30];

c = [0.000996388
0.000867473
0.000778297
0.000607295
0.00047941
0.000314246
0.000223062
6.98718E-05]; %SA
c(1:end,2) = 0.002916388 - (0.000996388-c); %AA

theta0 = [41.025*2; 0.000996388;	0.002916388]
options = optimset('FinDiffRelStep',1e-3,'DiffMinChange',1e-3,'DiffMaxChange',1);
[theta]=lsqcurvefit(@kinetics,theta0,t,c,[],[],options);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d',1:size(c,2)), 'Location','best')

end

.

[maurits]

Ah thanks very much! What was wrong with my code that it didn't iterate?

[Star Strider]

My pleasure!

‘What was wrong with my code that it didn't iterate?’

Most likely that you weren’t also estimating the initial conditions, since that’s the only significant change I made in your code (the other changes were simply asethetic). The initial concentrations should be considered parameters as well, since they’re likely measured (or at least should be measured rather than assumed), so they’re actually parameters and should be estimated as such.

If my Answer helped you solve your problem, please Accept it!

.