Function with multiple input parameters to be determined through fitting

Question

Jack 2024-6-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2126521-function-with-multiple-input-parameters-to-be-determined-through-fitting

评论： Jack 2024-6-8

Hi, I have a function to be fitted to some experimental data, but this function has multiple fitting parameters (6 parameters). I know that one of the ways to find these fitting parameters is to pass the function, the inital guess to the parameters and the lower and upper bounds into the lsqcurvefit function.

I have already created the function, but it seems like the software doesn't know that the p(1),p(2), p(3)... are parameters. Sorry if the function is very long.

What I wanted was p to be a vector with 6 elements p(1) p(2) p(3) p(4) p(5) p(6) but it seems like the software cant differentiate betwen p and p(1) p(2) p(3) p(4) p(5) p(6). When I try running the code, the error message I got was "Not enough input arguments"

Thank you for helping

function F = EM_SS(p, E_p)
F = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*(int(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))))), E, 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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jack 2024-6-7

please refer to this picture for the function to be fitted. everything other than

and π is a fitting parameter.

is given by:

where R is also a fitting parameter. There should be a total pf 6 fitting parameters.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-6-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2126521-function-with-multiple-input-parameters-to-be-determined-through-fitting#answer_1468996

编辑：Torsten 2024-6-8

在 MATLAB Online 中打开

You have to use "integral" instead of "int" and loop over the elements in E_p:

EM_SS([1 1 1 1 1 1],1)
ans = 3.2940e+03
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

15 个评论
显示 13更早的评论隐藏 13更早的评论

Jack 2024-6-8

编辑：Jack 2024-6-8

在 MATLAB Online 中打开

% Data files 
% global T_sample T_substrate Lambda E_p Abs
load Steady_State_Data.mat
h = 4.135667696*10^-15;
c = 3*10^8;
E_p = (Lambda*10^-9/h*c).^-1;
Abs = log10(T_substrate./T_sample);
lb = []; ub = [];
p0 = [0.13 0.1 1.6 0.05 3 1]; % initial guess for parameters (refer to line 13 for the order of parameters) 
% p0 = [A_1 A_2 E_g E_b R g]
% choose between different algorithm (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');
EM_SS([1 1 1 1 1 1], 1)
% solver
[p, fval, resnorm, residual, exitflag, jacobian] = lsqcurvefit(EM_SS, p0, E_p, Abs, lb, ub); 
plot(E_p, Abs) 

For reference, this is the script that runs the lsqcurefit function. and the EM_SS script is excatly similar to the one that you've editted for me.

T_sample, T_substrate and E_p are arrays of size 167 by 1

Torsten 2024-6-8

编辑：Torsten 2024-6-8

在 MATLAB Online 中打开

Steady_State_Data.mat

May I ask what is the difference between the 2? What effect does the capital letter have?

I use e_p as the array of hw and E_p as the i-th value from the array.

This code works, but the parameter vector becomes complex in the course of the computation which makes the "integral" function stop:

% Data files 
% global T_sample T_substrate Lambda E_p Abs
load Steady_State_Data.mat
h = 4.135667696*10^-15;
c = 3*10^8;
E_p = (Lambda*10^-9/h*c).^-1;
Abs = log10(T_substrate./T_sample);
lb = []; ub = [];
p0 = [0.13 0.1 1.6 0.05 3 1]; % initial guess for parameters (refer to line 13 for the order of parameters) 
% p0 = [A_1 A_2 E_g E_b R g]
% choose between different algorithm (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');
% solver
[p, fval, resnorm, residual, exitflag, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub); 
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 321)
        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 306)
    lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...
plot(E_p, Abs) 
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

Torsten 2024-6-8

I think the problems will be too problem-specific to find general ressources. As said: less parameters and good initial guesses for them would be helpful. Did you plot Abs against E_p for the initial guesses you provided to get an impression of the "initial fitting quality" ?

Jack 2024-6-8