lsqcurvefit "Function value and YDATA sizes are not equal" for a complicated fitting function
13 views (last 30 days)
Show older comments
Hi all, I have some experimental data to be fitted to a function that is quite complicated (refer to the attached picture). And I kept getting the message that the function size and YDATA size are not equal. How do I solve this?
clear; clc
load Steady_State_Data.mat % this contains the wavelength of light and absorbance of substrate and sample
% 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);
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
end
% Initial parameter guess and bounds
lb = []; ub = [];
p0 = [0.13 0.1 1.6 0.05 3 1]; % refer to the next line for their order
% p0 = [A1 A2 Eg Eb R g]
% choose between different algorithm of lsqcurvefit (3C1, comment those lines that are not choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt');
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective');
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point');
% fminunc
% optim_fminunc = optimsoptions('fminunc', 'Algorithm', 'Quasi-Newton');
% solver
[p, resnorm, residual, exitflag, output, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs);
% p = lsqcurvefit(@EM_SS, p0, E_p, Abs);
% Plot command
plot(E_p, Abs, 'o')
hold on
% plot(E_p, 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))))
plot(E_p, EM_SS(E_p, p))
xlabel('Probe energy (eV)')
ylabel('Absorbance.O.D')
legend('Experimental Data', 'Fitted Curve')
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);
6 Comments
Torsten
on 10 Jun 2024
Edited: Torsten
on 10 Jun 2024
Any idea how to resolve the issue of "Contour endpoints and waypoints must be finite."?
We already had this, too, in your previous question: the parameters and thus the lower limit of integration becomes complex-valued. Most probably, you somewhere take the square root of a negative number because "lsqcurvefit" suggests a negative parameter. In the call to "lsqcurvefit", use "lb" and "ub" to restrict the parameters to senseful values.
Answers (0)
See Also
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!