f =

Hello Yohay,

In the "fittype" there should be a parenthesis around "K_B*T" in the first term.

boltzmann = fittype...

( @(T,v) (m/(K_B*T))*v.*exp(((-m)/(2*K_B*T))*v.^2)...

, 'coefficient', 'T'...

, 'independent', 'v') ;

163 views (last 30 days)

Show older comments

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);

Torsten
on 18 Jul 2024 at 9:55

Edited: Torsten
on 18 Jul 2024 at 10:01

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)

int(f,x,0,Inf)

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)

sai charan sampara
on 18 Jul 2024 at 11:01

Edited: sai charan sampara
on 18 Jul 2024 at 11:02

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)

hist_y=hist_y/( trapz(hist_x,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]

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);

Alan Stevens
on 18 Jul 2024 at 7:56

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)

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

sai charan sampara
on 18 Jul 2024 at 9:17

Torsten
on 18 Jul 2024 at 14:11

Edited: Torsten
on 18 Jul 2024 at 14:13

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))

figure(1)

hold on

plot(hist_x,hist_y,'o');

plot(hist_x,f(sol,hist_x))

hold off

grid on

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

Start Hunting!