The parameters aren’t complex. All the results are real.
The problem with your code is that in MATLAB, function arguments and outputs are positional (except for the name-value pair arguments), and the output of your lsqcurvefit call did not account for all the outputs in order, so you weren’t getting what you thought you were getting.
I provided variable names for all of them here, for the ones you don’t need in your code, you can use the tilde (~) instead. An alternative to my changes would be:
[p, ~, residual, exitflag, output, ~, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub);
With that change, it now works —
clear; clc
load Steady_State_Data.mat % this contains the wavelength of light and absorbance of substrate and sample
% whos('-file','Steady_State_Data')
%% Preamble
% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
% Clean up of data to select range of values
Absorbance = log10(T_substrate./T_sample);
E = (h*c./(Lambda*10^-9));
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero
n = length(E) + 1;
% Data for fitting
E_p = E(n-N:167);
Abs = Absorbance(n-N:167) - min(Absorbance);
% Abs = Absorbance(n-N:167);
% Fitting function
function F = EM_SS(p, e_p)
for i = 1:numel(e_p)
E_p = e_p(i);
F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
(integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
p(2)*(2*pi*p(4)^3/2)*1/p(6)*(...
(1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
(1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
(1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
(1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
(1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
(1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
(1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6)));
end
F = F(:);
end
%% Curve fitting options
% Initial parameter guess and bounds
lb = [0,0,0,0,0,0]; ub = [Inf, Inf, 1.65, 0.03, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]
% lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',1000, 'MaxIterations', 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% Solver for lsqcurvefit
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub);
% p
% residual
% jacobian
%% Standard error calculation
ci = nlparci(p, residual, 'Jacobian', jacobian);
% lowerbar = p - ci(:,1);
% uppperbar = ci(:,2) - p;
Parameter_Confidence_Intervals = table(ci(:,1), p(:), ci(:,2),'VariableNames',{'CI 0.025','p','CI 0.975'})
[Ypred,delta] = nlpredci(@EM_SS,E_p,p,residual,'Jacobian',jacobian); % Prediction Confidence Limits
% if p == real(p)
% disp('Real')
% else
% disp('Complex with non-zero imaginary')
% end
%
% if residual == real(residual)
% disp('Real')
% else
% disp('Complex with non-zero imaginary')
% end
%
% if J == real(J)
% disp('Real')
% else
% disp('Complex with non-zero imaginary')
% end
%% Plot command
plot(E_p, Abs, '.')
hold on
plot(E_p, Ypred, '-r')
errorbar(E_p, Ypred, delta, 'r')
xlabel('Probe energy (eV)')
ylabel('Absorbance.O.D')
legend('Experimental Data', 'Fitted Curve', 'Prediction Confidence Intervals', 'Location','best')
hold off
%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
X1 = [' A1 = ', num2str(p1)];
X2 = [' A2 = ', num2str(p2)];
X3 = [' Eg = ', num2str(p3)];
X4 = [' Eb = ', num2str(p4)];
X5 = [' R = ', num2str(p5)];
X6 = [' g = ', num2str(p6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);
The second parameter is not significantly different from zero (the confidence limits are of opposite signs, so include zero), so it may not be necessary in your model. An alternative approach, if you have the Global Optimization Toolbox, would be to use the ga (genetic algorithm) function (my favourite, or one of the other multistart functions) first to search for the global minimum, and then use lsqcurvefit with those initial parameter estimates to refine the parameter estimates.
EDIT — (12 Jun 2024 at 05:47)
I just remembered that you want error bars, and the parameter confidence limits, while important, will not produce them. I added the nlpredci call to calculate the confidence intervals on the predictions (here, the fitted curve), and then used the errorbar function to plot them. (The error bars are vanishingly small in the plot, so you may need to zoom in to see them distinctly.)
.
