Something doesn't work for me in fit

5 次查看（过去 30 天）
Yohay
Yohay 2024-7-18
编辑： Torsten 2024-7-23
Hello friends.
A little introduction:
I do a simulation for particles moving inside a box, as part of the simulation I measure the speed of the particles, and the idea is that in the end I get the Maxwell-Boltzmann distribution of the speed.
I wanted to do a fit to the histogram of the velocity, according to the Boltzmann equation, so that I could find the temperature of the system, but for some reason the fit does not fit exactly on the data. it's similar, but not enough.
this is the fit that I get:
I'm uploading code with the histogram data, so you can run it yourself and see what you get.
I only change the histogram() in the figure, to plot(), so that it matches the data I uploaded.
This is the equation of Boltzmann distribution that I want to fit:
I will be happy if you can understand what is the problems..thanks!
This is the code, you can simple run it and see what is happening:
%my data:
hist_x=[10 30 50 70 90 110 130 150 170 190 210 230 250,...
270 290 310 330 350 370 390 410 430 450 470 490];
hist_y=[20 60 83 88 108 117 113 94 77 63 65 51 26,...
19 9 3 1 2 0 0 0 0 0 0 1];
%fit the velocity to Boltzmann distribution:
%parameters:
K_B=1.38e-23; %[m^2Kg/s^2K]
m=39.948*1.660e-27 ; %mass of Ar, convert from atomic mass to Kg [Kg]
boltzmann = fittype...
( @(T,v) (m/K_B*T)*v.*exp(((-m)/(2*K_B*T))*v.^2)...
, 'coefficient', 'T'...
, 'independent', 'v') ;
[fitted_curve] = fit(hist_x',hist_y',boltzmann, 'startPoint', 100) ;
%plot histogram and fit:
figure;
plot(hist_x,hist_y);
hold on
plot(fitted_curve)
title('histogram of velocity and fitted curve - Boltzmann distribution')
legend('histogram', 'fitted curve');
xlabel('velovity [m/s]')
ylabel('particles amount')
disp(fitted_curve);
  9 个评论
显示 7更早的评论
Torsten
Torsten 2024-7-18
编辑：Torsten 2024-7-18
Ran in:
@sai charan sampara
Your function is not a distribution - it doesn't integrate to 1.
syms a x real positive
f = 1/sym(pi)^0.5*a^0.5*x*exp(-a*x^2)
f = 
int(f,x,0,Inf)
ans = 
@Yohay
Your data don't integrate to 1. So it's not possible to get a good fit.
hist_x=[10 30 50 70 90 110 130 150 170 190 210 230 250,...
270 290 310 330 350 370 390 410 430 450 470 490];
hist_y=[20 60 83 88 108 117 113 94 77 63 65 51 26,...
19 9 3 1 2 0 0 0 0 0 0 1];
hist_y=hist_y/sum(hist_y);
trapz(hist_x,hist_y)
ans = 19.7900
sai charan sampara
sai charan sampara 2024-7-18
编辑：sai charan sampara 2024-7-18
Ran in:
You are right Torsten. I didn't think of that. Will scaling the data by the above integral value to make the area become 1 and then doing a fit be a valid approach?
hist_x=[10 30 50 70 90 110 130 150 170 190 210 230 250,...
270 290 310 330 350 370 390 410 430 450 470 490];
hist_y=[20 60 83 88 108 117 113 94 77 63 65 51 26,...
19 9 3 1 2 0 0 0 0 0 0 1];
hist_y=hist_y/(sum(hist_y));
trapz(hist_x,hist_y)
ans = 19.7900
hist_y=hist_y/( trapz(hist_x,hist_y));
trapz(hist_x,hist_y)
ans = 1
K_B=1.38e-23; %[m^2Kg/s^2K]
m=39.948*1.660e-27 ; %mass of Ar, convert from atomic mass to Kg [Kg]
boltzmann = fittype...
( @(T,v) (m/(K_B*T))*v.*exp(((-m)/(2*K_B*T))*v.^2)...
, 'coefficient', 'T'...
, 'independent', 'v') ;
[fitted_curve] = fit(hist_x',hist_y',boltzmann, 'startPoint', 100) ;
%plot histogram and fit:
figure;
plot(hist_x,hist_y);
hold on
plot(fitted_curve)
disp(fitted_curve);
General model: fitted_curve(v) = (m/(K_B*T))*v.*exp(((-m)/(2*K_B*T))*v.^2) Coefficients (with 95% confidence bounds): T = 100 (71.99, 128)

请先登录，再进行评论。

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

回答（2 个）

Alan Stevens
Alan Stevens 2024-7-18
Ran in:
Here's a fit using fminsearch (I just used a quick but coarse approach- the code can be made much cleaner!).
v=[10 30 50 70 90 110 130 150 170 190 210 230 250,...
270 290 310 330 350 370 390 410 430 450 470 490];
y=[20 60 83 88 108 117 113 94 77 63 65 51 26,...
19 9 3 1 2 0 0 0 0 0 0 1];
bar(v,y,'c')
hold on
X0 = [1 -0.01];
X = fminsearch(@fn, X0);
vv = linspace(0,500);
Y = X(1)*vv.*exp(X(2)*vv.^2);
plot(vv,Y,'k-')
xlabel('v'), ylabel('y')
format long
disp(X)
1.742146355175938 -0.000043540164592
function Z = fn(X)
v=[10 30 50 70 90 110 130 150 170 190 210 230 250,...
270 290 310 330 350 370 390 410 430 450 470 490];
y=[20 60 83 88 108 117 113 94 77 63 65 51 26,...
19 9 3 1 2 0 0 0 0 0 0 1];
MB = X(1)*v.*exp(X(2)*v.^2);
Z = norm(MB-y);
end
  1 个评论
sai charan sampara
sai charan sampara 2024-7-18
There is only one variable/coefficient "T" that is varying. "X(1)" and "X(2)" are dependent quantities.

请先登录，再进行评论。

Torsten
Torsten 2024-7-18
编辑：Torsten 2024-7-18
Ran in:
hist_x=[10 30 50 70 90 110 130 150 170 190 210 230 250,...
270 290 310 330 350 370 390 410 430 450 470 490];
hist_y=[20 60 83 88 108 117 113 94 77 63 65 51 26,...
19 9 3 1 2 0 0 0 0 0 0 1];
hist_y=hist_y/trapz(hist_x,hist_y);
K_B=1.38e-23; %[m^2Kg/s^2K]
m=39.948*1.660e-27 ; %mass of Ar, convert from atomic mass to Kg [Kg]
f = @(T,v)m./(K_B*T).*v.*exp(-m./(2*K_B*T)*v.^2);
T0 = 100;
sol = lsqcurvefit(f,T0,hist_x,hist_y,[],[],optimset('TolX',1e-10,'TolFun',1e-10))
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
sol = 54.8456
figure(1)
hold on
plot(hist_x,hist_y,'o');
plot(hist_x,f(sol,hist_x))
hold off
grid on
  3 个评论
显示 1更早的评论
Yohay
Yohay 2024-7-23
@Torsten
Hi, thanks!
How you get this parametrs? you do this on the curve fit app?
Torsten
Torsten 2024-7-23
编辑：Torsten 2024-7-23
Copy the code in your MATLAB editor and execute it.
As you can see from the output above, T comes out to be approximately 55.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by