lsqcurvefit (with system of ODEs) does not iterate, it only displays the initial guess value
12 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Star Strider
on 26 Jun 2023
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
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
.
2 Comments
Star Strider
on 27 Jun 2023
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!
.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!