How can I get this code to work? I am trying to fit the equation that contains step and delta function to experimental data. But I am getting the following error.

3 ビュー (過去 30 日間)
Bibek Dhami
Bibek Dhami 2022 年 6 月 24 日
コメント済み: Bibek Dhami 2022 年 6 月 30 日
A=xlsread('Data');
xdata= A(:,1);
ydata = A(:,2) ;
plot(xdata,ydata,'r.','markersize',30);
set(gca,'XDir','reverse');
box on
options =optimoptions(@lsqcurvefit,'Algorithm','trust-region-reflective','StepTolerance',1e-19,'MaxFunctionEvaluations',1e10)
options =
lsqcurvefit options: Options used by current Algorithm ('trust-region-reflective'): (Other available algorithms: 'levenberg-marquardt') Set properties: Algorithm: 'trust-region-reflective' MaxFunctionEvaluations: 1.0000e+10 StepTolerance: 1.0000e-19 Default properties: CheckGradients: 0 Display: 'final' FiniteDifferenceStepSize: 'sqrt(eps)' FiniteDifferenceType: 'forward' FunctionTolerance: 1.0000e-06 JacobianMultiplyFcn: [] MaxIterations: 400 OptimalityTolerance: 1.0000e-06 OutputFcn: [] PlotFcn: [] SpecifyObjectiveGradient: 0 SubproblemAlgorithm: 'factorization' TypicalX: 'ones(numberOfVariables,1)' UseParallel: 0
MI=lsqcurvefit(@model,[0.2,0.0403,2.03],xdata,ydata,[],[],options);
Error using lsqncommon
Objective function is returning undefined values at initial point. lsqcurvefit cannot continue.

Error in lsqcurvefit (line 295)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...
fitdata=model(MI,xdata);
hold on
plot(xdata, fitdata,'k--','linewidth',2);
function alpha=model(var,energy)
alpha = var(1)*heaviside(energy-var(3)).*((pi.*exp(pi.*energy))./sinh(pi.*energy))+...
dirac(energy-var(3)+(var(2)))+dirac(energy-var(3)+(var(2)./4))+...
dirac(energy-var(3)+(var(2)./9));
end
  1 件のコメント
dpb
dpb 2022 年 6 月 24 日
>> F=model([0.2,0.043,2.03],A(:,1));
'heaviside' requires Symbolic Math Toolbox.
Error in model (line 2)
alpha = x(1)*heaviside(energy-x(3)).*((pi.*exp(pi.*energy))./sinh(pi.*energy))+...
>>
so can't run your function here...first thing to do, however, is to debug your model function before trying it inside lsqcurvefit
The error is telling you something ain't right there -- you've not returning values that are defined and finite at the start point...

サインインしてコメントする。

採用された回答

Walter Roberson
Walter Roberson 2022 年 6 月 24 日
Your Energy are in the 500 range. exp(pi*energy) overflows double precision. Your code then divides by large values, but you have already overflowed.
Ideally your code should be revised to deliberately avoid overflowing. For example, instead of exp(pi*energy)/sinh(something) you could try exp(pi*energy - log(sinh(something)) expecting the sinh to be large enough for the log to be comparable to pi*energy.
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/1044920/Data.xlsx';
A = readmatrix(filename);
xdata = A(:,1);
ydata = A(:,2) ;
plot(xdata, ydata, 'r.', 'markersize', 30);
set(gca,'XDir','reverse');
box on
options =optimoptions(@lsqcurvefit,'Algorithm','trust-region-reflective','StepTolerance',1e-19,'MaxFunctionEvaluations',1e10)
options =
lsqcurvefit options: Options used by current Algorithm ('trust-region-reflective'): (Other available algorithms: 'levenberg-marquardt') Set properties: Algorithm: 'trust-region-reflective' MaxFunctionEvaluations: 1.0000e+10 StepTolerance: 1.0000e-19 Default properties: CheckGradients: 0 Display: 'final' FiniteDifferenceStepSize: 'sqrt(eps)' FiniteDifferenceType: 'forward' FunctionTolerance: 1.0000e-06 JacobianMultiplyFcn: [] MaxIterations: 400 OptimalityTolerance: 1.0000e-06 OutputFcn: [] PlotFcn: [] SpecifyObjectiveGradient: 0 SubproblemAlgorithm: 'factorization' TypicalX: 'ones(numberOfVariables,1)' UseParallel: 0
MI = lsqcurvefit(@model,[0.2,0.0403,2.03],xdata,ydata,[],[],options);
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
fitdata = model(MI,xdata);
hold on
plot(xdata, fitdata,'k--','linewidth',2);
function alpha = model(var,Energy)
energy = sym(Energy);
alpha = var(1)*heaviside(energy-var(3)).*((pi.*exp(pi.*energy))./sinh(pi.*energy))+...
dirac(energy-var(3)+(var(2)))+dirac(energy-var(3)+(var(2)./4))+...
dirac(energy-var(3)+(var(2)./9));
alpha = double(alpha);
end
  4 件のコメント

サインインしてコメントする。

その他の回答 (1 件)

Cris LaPierre
Cris LaPierre 2022 年 6 月 24 日
編集済み: Cris LaPierre 2022 年 6 月 27 日
The problem is that your objective function is returning NaN values. It is specifically the exp and sinh parts of the equation. You may need to consider the scale of your x values.
A=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1044920/Data.xlsx')
A = 83×2
496.7700 0.2983 497.3600 0.3010 497.9500 0.3000 498.5400 0.3009 499.1300 0.3037 499.7200 0.3059 500.3000 0.3065 500.8900 0.3088 501.4800 0.3101 502.0700 0.3113
xdata= A(:,1);
ydata = A(:,2);
var = [0.2,0.0403,2.03];
energy = xdata;
% exp part retuning inf
pi.*exp(pi.*energy)
ans = 83×1
Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
% sinh part returning inf
sinh(pi.*energy)
ans = 83×1
Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf
% dividing an inf by inf returns a nan
((pi.*exp(pi.*energy))./sinh(pi.*energy))
ans = 83×1
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  5 件のコメント
Bibek Dhami
Bibek Dhami 2022 年 6 月 27 日
HI @Walter Roberson I just used your technique but the code keeps running without giving output. I need to stop the code after 5 minutes.

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by