lsqcurvefit "Function value and YDATA sizes are not equal" for a complicated fitting function

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);

Error using lsqcurvefit (line 312)
Function value and YDATA sizes are not equal.

% 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);

[Torsten]

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]

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')

  Name               Size            Bytes  Class     Attributes

  Absorbance       167x1              1336  double              
  E                167x1              1336  double              
  Lambda           167x1              1336  double              
  T_sample         167x1              1336  double              
  T_substrate      167x1              1336  double              

% 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);

Error using integralCalc (line 33)
Contour endpoints and waypoints must be finite.

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in solution>EM_SS (line 26)
(integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...

Error in lsqcurvefit/objective (line 337)
F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});

Error in snls (line 344)
newfvec = feval(funfcn{3},xcurr,varargin{:});

Error in lsqncommon (line 204)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,F,Jac,caller, ...

Error in lsqcurvefit (line 322)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...

% 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]

@Torsten adding the line F = F(:) will make the function loop infinitely. I've actually tried without the F = F(:) and it worked. But just as I change my dataset, the code now doesn't work anymore

[Jack]

Any idea how to resolve the issue of "Contour endpoints and waypoints must be finite."? Should I change the limits of the integration to a finite value instead of infinity?

[Torsten]

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]

I see I see. Thank you all for the help. Appreciate it.