kill(all)$

Maxwell - Boltzmann distribution .(appendix 2),version 1.0
peter.vlasschaert@gmail.com,23/02/2019

see : website
Maxwell - Boltzmann distribution .(appendix 1),version 1.0

load packages
=============

load(lrats)$

find : equilibrium distribution ( = most general probability function).

Def  : n_i = g_i *exp(-α-β*W_i)

i       = ith energy level
g_i  = probability different for every energy level  ( see :appendix 1)
α, β = lagrange multiplicators (from the constraint ,see appendix 1)

H ow  : Max (P)  better do calculation  Max(log(P))
=========================================

1e)  (stir)ling approximation : lnx! ( Maxima  : log (x!))
----------------------------------------------------------------------------
because : number n_i very big numbers ( = number of particles)

p1:  log(x!) = integrate(log(x),x);

p2 : subst([m=k,x=m],p1);

construct function : 'p1'
====================

p2a:stir(m) := (l:subst([m=k,x=m],p1));

stir(vv);

2e)  example : three energy level (Max(P) => Max(log(P))).
-------------------------------------------------------------------------------------
use : ' p7d ' appendix : 1

p3:P=product(g[i]^n[i],i,1,3)/(n[1]!*n[2]!*n[3]!);

p3a:log(p3);

p3b :lhs(p3a)=rhs(p3a), logexpand=super;

3e)  from constraint  'p1' and 'p2'  appendix 1
--------------------------------------------------------------------
example : three energy levels  (k =3 ,appendix 1)

p4:     N =sum(n[i],i,1,3);
p4a:  W= sum(n[i]*W[i],i,1,3);

How to solve problem.
===================

a)   n_i  : pos integer
------------------------------

assume(n[1]>0,n[2]>0,n[3]>0);

b)  approximation from stirling (use 'p3b' above)
-------------------------------------------------------------------------

stir(n[1]);

recursive :
=========

p5:for i:1 thru 3 do(
p3b:subst([stir(n[i])],p3b))$

p5a:p3b;

p5b :logcontract(p5a);

p5c :expand(p5b);

c)  diff(log(P))
---------------------

p6:diff(p5c);

recursive :   del(g[i]) = 0 , i=1,2,3
=========

p6a:for i:1 thru 3 do(
p6:subst([del(g[i]) =0],p6))$

p6b:p6;

?  :del(P) = 0

p6c:lhs(p6b) = -'sum(log(n[i]/g[i])*del(n[i]),i,1,3);

p6b:p6c*(-1);

p6b1:subst([del(P) = 0],p6b);

rem : del(n_i) = Δ n_i  different form zero (change)

lhs :
====

p6c:subst([del(P) = 0],lhs(p6b));

rhs :
====

p6d1:part(rhs(p6b),0);
p6d2:part(rhs(p6b),1);
p6d3:part(p6d2,1);
p6d4:part(p6d2,2);

result :  'del(n_i)' ≠ 0
==============

p7: p6d3 = p6c;

constraint :  α,β
=============

a) constraint : 'p1', appendix 1
----------------------------------------------

p8:   dN = 'sum(del(n[i]),i,1,3);

p8a:α*(p8);

because :  N = constant  , dN = 0

p8b:subst([dN = 0],p8a);

p8c:intosum(p8b);

b) constraint : 'p2', appendix 1
----------------------------------------------

p9:0='sum(W[i]*del(n[i]),i,1,3);
p9a:β*(p9);

p9b:intosum(p9a);

L = function + α (first constraint) + β(second constraint)
==============================================

p10:p6b1+p8c+p9b;

p10a:lhs(p10) =part(rhs(p6b1),1) + part(rhs(p8c),1) + part(rhs(p9b),1);

p10b:lhs(p10a) = coeff(rhs(p10a),del(n[i]))*del(n[i]);

total result : 'Def'
==============

p10c: lhs(p10b)  = part(rhs(p10b),1);

p10b:solve(p10c,n[i])[1];

i th : level  for the energy.

'-----------------------------------------------------------------------------------------------------'