Something doesn't work for me in fit

4 views (last 30 days)
Yohay
Yohay on 18 Jul 2024
Edited: Torsten on 23 Jul 2024
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 Comments
Torsten
Torsten on 18 Jul 2024
Edited: Torsten on 18 Jul 2024
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 = 
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 on 18 Jul 2024
Edited: sai charan sampara on 18 Jul 2024
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)

Sign in to comment.

Answers (2)

Alan Stevens
Alan Stevens on 18 Jul 2024
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 Comment
sai charan sampara
sai charan sampara on 18 Jul 2024
There is only one variable/coefficient "T" that is varying. "X(1)" and "X(2)" are dependent quantities.

Sign in to comment.


Torsten
Torsten on 18 Jul 2024
Edited: Torsten on 18 Jul 2024
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 Comments
Yohay
Yohay on 23 Jul 2024
Hi, thanks!
How you get this parametrs? you do this on the curve fit app?
Torsten
Torsten on 23 Jul 2024
Edited: Torsten on 23 Jul 2024
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.

Sign in to comment.

Categories

Find more on Linear and Nonlinear Regression in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!