boltzmann shot 1

An Attempt at the Boltzmann Problem PhL 6.2..15 Here is the problem as stated in my blackbody doc (but perhaps this is not really the blackbody problem) Ω(N1,N2...Nm) = .... = number of "accessible states" for the full set of particles Σi=1m Ni = M // number of particles is a constant (so not photons I guess) Σi=1m Ni εi = U // total energy is a constant (adiabatic I guess) Here there are a finite set m of discrete energy states called εi. Each state has Ni particles in it. Each state has degeneracy gi which only effects Ω. You can I guess have Ni larger than gi if you want, though I think typically Ni < gi. The total number of particles is M, so you cannot create or destroy particles. The total energy of the particles is U. So this is a classic Lagrange Multiplier problem! I will try to solve it using my new LM document. Let's restate the problem using Ni → xi to fit into my format f(x1,x2...xm) = .... Σi=1m xi - M = 0 a(x1,x2...xm) = 0 Σi=1m xi εi - U = 0 b(x1,x2...xm) = 0 Well, right away we have a distinction with the classical LM case. Here the xi are all positive integers, not real numbers! But I guess I can "fill in" using the gamma function and x! = Γ(x+1). Then ∂x (x!) = ∂x Γ(x+1) = [ ∂yΓ(y)]y=x+1 = ? From GR p 902 we know that This says ψ(x) = ∂xln(Γ(x)) = 1/Γ * ∂xΓ(x) ∂xΓ(x) = Γ(x)ψ(x) so this gets us involved with "the psi function" which is also known as the "digamma function". Start with Ω(N1,N2...Nm) = .... lnΩ = ( N1lng1 - ln N1! ) + ( N2lng2 - ln N2! ) + .... + ( Nmlngm - ln Nm! ) = Σi=1m [ Ni ln gi ] - Σi=1m [ ln Ni! ] Perhaps lnΩ is an easier function to maximize. Stirling says ln x! = x ln x - x. , so let's assume that all the Ni are very large numbers and then ln Ni! ≈Ni ln Ni - Ni. Then we find f ≡ lnΩ ≈ Σi=1m [ Ni ln gi ] - Σi=1m [ Ni ln Ni - Ni ] = Σi=1m [ Ni ln ] + Σi=1m Ni = Σi=1m [ Ni ln ] + M // Zemansky (10-6) and this certainly is a simpler functional form. Now when we compute the partials fi only one term is picked out, so we have fi = ∂i [ Ni ln ] = ? Simplify for the moment and write ∂x [ x ln ] = ∂x [ x ln g - x ln x ] = ln g - ∂x [x ln x] = ln g - { x * 1/x + ln x } = ln g - {1 + ln x } = ln (g/x) - 1 which is a very simple result! Then fi = ∂i [ Ni ln ] = ln - 1 This may be wrong because I have used the constraint Σi=1m Ni = M, but fi should be computed without applying any of the constraints. Thus, the correct fi will be different from what I show here. Let's go start a new doc boltz shot 2. So now we restate our LM problem assuming we are only interested in regions where all Ni are very large: f(N1,N2...Nm) = Σi=1m [ Ni ln ] + M // maximize f a(N1,N2...Nm) = Σi=1m Ni - M = 0 b(N1,N2...Nm) = Σi=1m Ni εi - U = 0 or in our usual variable names f(x1,x2...xm) = Σi=1m [ xi ln ] + M // maximize f a(x1,x2...xm) = Σi=1m xi - M = 0 b(x1,x2...xm) = Σi=1m xi εi - U = 0 We know at once that fi = ln - 1 = ln(gi/xi) - 1 ai = 1 bi = εi That seems to make a relatively simple R matrix: ln(g1/x1) - 1 ln(g2/x2) - 1 ..... ln(gm/xm) - 1 R = 1 1 ..... 1 ε1 ε2 ..... εm Now follow the algorithm and see if it leads anywhere useful: H(r,λ) ≡ f(r) + λ1a(r) + λ2 b(r) + ...... + λS-1 q(r) . (2.1) or H(r,λ) ≡ Σi=1m [ xi ln ] + M + λ1[ Σi=1m xi - M] + λ2[ Σi=1m xi εi - U] Next,. Hi = fi + λ1ai + λ2 bi+ ...... + λS-1 qi i = 1,2....N (2.2) becomes Hi = ln(gi/xi) - 1 + λ1 + λ2 εi Now set all partials to zero, so we then have ln(gi/xi) - 1 + λ1 + λ2 εi = 0 i = 1,2....m Here then is the set of equations we have to solve: ln(gi/xi) - 1 + λ1 + λ2 εi = 0 i = 1,2....m Σi=1m xi = M Σi=1m xi εi = U There are m+2 equations in m+2 unknowns. Let's try to find candidate solutions for the λi as described in my LM doc: [ not the right path to take here, but here it is anyway ] - = (2.7) - = (2.7) Now I pick various choices of 2x2 from this. For example = or = = -1 Maple computes the inverse matrix: Thus, M-1 = and then = and this says λ1 = [ -εjfi + εifj] λ2 = [ fi - fj] Not bad. Now fi - fj = [ln(gi/xi) - 1] - [ln(gj/xj) - 1] = [ ln(gi/xi) - ln(gj/xj) ] And we can write all out as λ1 = { -εj [ln(gi/xi) - 1] + εi[ln(gj/xj) - 1]} = { (εj-εi) - εjln(gi/xi) + εiln(gj/xj)} = -1 + [εiln(gj/xj) - εjln(gi/xi)] λ2 = [ ln(gi/xi) - ln(gj/xj) ] I have no idea whether this will lead anywhere, but it is a great test of my LM doc. Now back to the equation set. I will mark the choice i,j values as I,J, so ln(gi/xi) - 1 + λ1 + λ2 εi = 0 i = 1,2....m Σi=1m xi = M Σi=1m xi εi = U becomes ln(gi/xi) - 2 + {εIln(gJ/xJ) - εJln(gJ/xI)} + [ ln(gI/xI) - ln(gJ/xJ) ] εi = 0 i = 1,2,....m Σi=1m xi = M Σi=1m xi εi = U So I have eliminated the λ1 and λ2. ALL candidate solutions must solve these equations where you make a selection of I≠J from 1...m. Now let's back up a little to ln(gi/xi) - 1 + λ1 + λ2 εi = 0 i = 1,2....m Σi=1m xi = M Σi=1m xi εi = U Remember that the εi are givens, not unknowns. So try to solve the first equation for xi ln(gi) - ln(xi) - 1 + λ1 + λ2 εi = 0 ln(xi) = ln(gi) - 1 + λ1 + λ2 εi ln(xi/gi) = - 1 + λ1 + λ2 εi xi = gi exp( - 1 + λ1 + λ2 εi) = gi e-1 e-λe-λε Once again, I have just solved the first m equations to find xi = gi e-(λ+1)e-λε and for the very first time I am seeing something looking like a "Boltzmann energy exponential"! OK, then we are left with these 2 equations in the two unknowns λ1 and λ2 : Σi=1m gi e-(λ+1)e-λε = M Σi=1m giεi e-(λ+1)e-λε = U or Σi=1m gie-λε = M e(λ+1) Σi=1m εigie-λε = Ue(λ+1) Divide these equations to get [Σi=1m εigie-λε] / [Σi=1m gie-λε] = U/M = u = average energy per particle or = u VERY interesting! Let's now define pi ≡ K gie-λε = probability that a particle is in state i. K = not yet known Then the above states = u and this makes total sense with that interpretation. This sort of came out of nowhere! Now solve for λ2 [Σi=1m εigie-λε] / [Σi=1m gie-λε] = U/M = u Σi=1m εigie-λε = u Σi=1m gie-λε Σi=1m (εi- u)gie-λε = 0 Here the εi are known, u is U/M which is also known, gi is known. So how solve for λ2? Σi ai exp(-xbi) = 0 solve for x Define z = e-x so this then says F(z) = Σi aizb = 0 What kind of equation is this? A polynomial with real instead of integer coefficients. I see the term "exponential polynomials" bandied about and I found one downloadable paper. I think for general bi this is a transcendental equation. If the bi are rational, you have a polynomial equation. Go back to Σi ai exp(-xbi) = 0 solve for x Write bi = Ni/Mi and find some large J = lcm(Mi) so that J = KiMi. Then exp(-xbi) = exp(-xNi/Mi) = exp(-xKiNi/J) = exp(-KiNi [x/J] ) = (e-x/J)KN Now let z ≡ e-x/J so you then have Σi ai zKN = 0 and then you have converted your exponential polynomial equation into a normal polynomial equation. In this case, the solution z is an algebraic function according to the wiki definition, since z is then the solution of a polynomial equation with coefficients which are polynomials (here just constants ai). If the bi are not rational, like b2 = π, then it is a transcendental equation. In any event, generally higher order polynomials have no analytic solutions for roots. at least no simple ones, and an example is given in this clip: So the upshot of this discussion is that we wanted to find λ2 from Σi=1m (εi- u)gie-λε = 0 and in general this needs a numeric solution. For completely general εi values, it is a transcendental equation, while for rational εi values it is a polynomial equation. So I won't be able to analytically solve for the value λ2. There might not even BE a solution, but I think in my case there is. Next question: Assuming we have a solution for λ2, can we solve for λ1 ? Our equations were (1/M) Σi=1m gie-λε = e(λ+1) = (1/M) Σipi (1/U) Σi=1m εigie-λε = e(λ+1) = (1/U) Σipiεi The first equation says Σipi = M e(λ+1) and this then gives a solution value for λ1 if you know λ2. The second equation they just says = u and so we can in fact satisfy both. Question: How does this relate do the supposed solution values of λ1 and λ2 I found using my LM document method, as outlined above: λ1 = -1 + [εiln(gj/xj) - εjln(gi/xi)] λ2 = [ ln(gi/xi) - ln(gj/xj) ] These are candidate solution pairs. But they are expressed in terms of xi and xj which are unknowns of the problem, so I guess there is no conflict. OK, I want now to continue on this general path to see where I have gotten! Question: How do I know I have found a maximum? f(x1,x2...xm) = Σi=1m [ xi ln ] + M I have NO IDEA how to show this is a max. I know that the solution I found is ln(xi/gi) = - 1 + λ1 + λ2 εi I guess I could install that to get fextrememum(x1,x2...xm) = Σi=1m [ - xi ln(xi/gi) ] + M = - Σi xi ln(xi/gi) + M = - Σi xi [- 1 + λ1 + λ2 εi] + M = - Σi xi [- 1 + λ1 + λ2 εi] + Σixi = - Σi xi [- 2 + λ1 + λ2 εi] = Σi xi [ 2 - λ1 - λ2 εi] = (2-λ1)Σi xi - λ2 Σi xiεi = (2-λ1)M - λ2U Well, maybe this is what you do. Somehow go to a small neighborhood of your solution, construct a legal dr, and see if f dr is positive or negative! A legal dr has to respect the two constraints, Σixi = M Σxi εi = U If m = 3, we have x1+ x2+ x3 = M and this is a plane in E3. Well OK, go back to a(x1,x2...xm) = Σixi - M = 0 a = (1,1...1) m components b(x1,x2...xm) = Σixi εi - U = 0 b = (ε1, ε2....εm) m components There are the normals to the constraint planes. Any legal dr must be perp to these two normals. So need dr a = 0 or Σdxi = 0 dM = 0 dr b = 0 or Σεidxi = 0 dU = 0 Let's try a Taylor expansion around the solution point f(x1,x2...xm) = Σi [ xi ln ] + M ∂f/∂xi = fi = ln - 1 // first derivatives ∂2f/∂xi2 = - 1/xi ∂2f/∂xi∂xj = -(1/xi) δi,j The surface is "cupping down" in all directions. This is just like my hemisphere example in LM doc, it just happens. That is encouraging. Is that all I really need to know?