lsqcurvefit "Function value and YDATA sizes are not equal" for a complicated fitting function
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
Didn't I answer this already in a previous question of yours ?
Add the line
F = F(:);
as last line in "EM_SS" to make F a column vector.
Star Strider
on 10 Jun 2024
Moved: Star Strider
on 10 Jun 2024
Since ‘E_p’ is a coilumn vector, ‘F’ must also be a column vector. That is easy to fix, however once done, that brings up a different problem in the integral call about non-finite contour points and waypoints. I defer to you for that solution.
clear; clc
load Steady_State_Data.mat % this contains the wavelength of light and absorbance of substrate and sample
whos('-file', 'Steady_State_Data')
% 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
F = F(:);
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);
.
Jack
on 10 Jun 2024
Jack
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.
Jack
on 11 Jun 2024
Answers (0)
Categories
Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange
on 11 Jun 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
