Fitting a blackbody curve with two variables - having issues with lsqcurvefit stopping and unsure how to use options effectively

11 次查看（过去 30 天）
Caelan
Caelan 2025-6-3
评论： Caelan 2025-6-10
Ran in:
Hello,
I am trying to fit a blackbody curve to some data I have and need to find the value for two variables. I have the following code:
format long g
%fit BlackBody to data
%constants needed for BB equation
% SI UNITS
c=2.998e8;
h=6.6261e-34;
kb=1.38e-23;
%wavelgnth (xdata) and flux densities (ydata).
xdata = [0.000000212, 0.000000231, 0.000000291, 3.465E-07, 3.543E-07, 0.000000365, 4.392E-07, 0.000000445, 0.000000477, 5.468E-07, 0.000000551, 6.231E-07, 0.000000658, 7.625E-07, 0.000000806, 9.134E-07, 0.00000122, 0.00000163, 0.00000215];
ydata = [1.32933E-13, 1.11645E-13, 1.04652E-13, 6.34698E-14, 5.64255E-14, 5.39367E-14, 3.80465E-14, 3.66716E-14, 3.55176E-14, 2.93087E-14, 3.10419E-14, 2.40891E-14, 2.09837E-14, 1.42909E-14, 1.2889E-14, 1.01056E-14, 5.81715E-15, 3.5086E-15, 1.99137E-15];
%plot figure with data
figure
plot(xdata, ydata, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 5);
xlabel('Wavelength (m)');
ylabel('Flux Density (wavelength depedence, W/m^3)');
%BB function - using variable names from MATLAB documentation for ease of
%learning to use fitting function. Here x(1) is the aplha multiplication
%factor and x(2) is temperature. xdata is wavelength
fun = @(x,xdata) x(1)*(2.*pi.*h.*(c^2)./(xdata.^5)).*(1./(exp((h.*c)./(xdata.*kb.*x(2)))-1));
%fitting
x0 = [1e-23 28000]; %starting point for fitting
x = lsqcurvefit(fun,x0,xdata,ydata)
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
x = 1×2
1.0e+00 * 1.87324626094452e-30 28000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
%plot fit
hold on
plot(xdata,fun(x,xdata),'r-');
The values produced for x show a value I would expect for x(2) but not for x(1). I would expect it to be around e-23. The solver is stopping, I have tried to edit the options for step count, but I'm unsure of what values to change within options and what would be reasonable numbers to use.
Anyone have any idea how to get a more accurate x(1) value and how to get the solver to run completely?
Thank you!
  1 个评论
Sam Chak
Sam Chak 2025-6-4
Hi @Caelan
What estimation formula did you use to determine that the value for the parameter x(1) is expected to be on the order of ? Are all blackbody radiations characterized in this manner?

请先登录，再进行评论。

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

采纳的回答

John D'Errico
John D'Errico 2025-6-3
Your data does not seem to fit terribly well to the model you propose.
c=2.998e8;
h=6.6261e-34;
kb=1.38e-23;
%wavelgnth (xdata) and flux densities (ydata).
xdata = [0.000000212, 0.000000231, 0.000000291, 3.465E-07, 3.543E-07, 0.000000365, 4.392E-07, 0.000000445, 0.000000477, 5.468E-07, 0.000000551, 6.231E-07, 0.000000658, 7.625E-07, 0.000000806, 9.134E-07, 0.00000122, 0.00000163, 0.00000215];
ydata = [1.32933E-13, 1.11645E-13, 1.04652E-13, 6.34698E-14, 5.64255E-14, 5.39367E-14, 3.80465E-14, 3.66716E-14, 3.55176E-14, 2.93087E-14, 3.10419E-14, 2.40891E-14, 2.09837E-14, 1.42909E-14, 1.2889E-14, 1.01056E-14, 5.81715E-15, 3.5086E-15, 1.99137E-15];
I'll use my own fminspleas to do the fitting, as it is far more robust to problems.
fun = @(x,xdata) (2.*pi.*h.*(c^2)./(xdata.^5)).*(1./(exp((h.*c)./(xdata.*kb.*x))-1));
[INLP,ILP] = fminspleas({fun},28000,xdata,ydata)
INLP =
14620.5497503281
ILP =
1.48319389723619e-29
plot(xdata,ydata,'o')
hold on
plot(xdata,ILP*fun(INLP,xdata),'-s')
Sigh. The plot did not come out terribly clearly. But you can see it chose a black body temperature around half of what you expected.
  2 个评论
Caelan
Caelan 2025-6-10
This was really helpful thanks! You're right, I had an error in the data, I have rectified this now and the code worked perfectly.
Caelan
Caelan 2025-6-10
Also - since this is your own function - can I ask if there is a way to determine the errors on the fitted INLP and ILP values?

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by