Re: FindRoot and Bose-Einstein distribution
- To: mathgroup at smc.vnet.net
- Subject: [mg82887] Re: FindRoot and Bose-Einstein distribution
- From: danl at wolfram.com
- Date: Fri, 2 Nov 2007 03:32:11 -0500 (EST)
- References: <fg6ql3$dmc$1@smc.vnet.net>
On Oct 30, 3:40 am, P_ter <peter_van_summe... at yahoo.co.uk> wrote:
> Hello,
> I have two functions and two values:
> Clear[p, q, k, mm, mn]
> n[p_, q_] := Sum[k/(Exp[p + k q ] - 1), {k, 1, 10000}]
> m[p_, q_] := Sum[1/(Exp[p + k q ] - 1), {k, 1, Infinity}]
> mm = 0.501
> mn = 41.0959
> The first two equations are the Bose-Einstein distribution. Given mm and mn, find p and q.
> A first estimation is: p=5.086,q= 0.01226
> I know that n[p_,q_] is stable until k = 50000 for p=5.086 and q= 0.01226. My check is: n[5.086,0.01226]= 41.1969
> m[p_,q_] is also stable until 100000 and m[5.086,0.01226]= 0.502762
It is not clear to me what is meant by "stable until k=...". Do you
actually intend to do infinite summations for both m and n, and are
checking to see how far finite truncations can be computed wo some
precision?
> Everything seems ok.
> So, I tried:
> FindRoot[{m[p, q] - mm, n[p, q] - mn}, {{p, 5}, {q, 0.1}}]
> The answer from FindRoot is that the Sum does not converge.
> What did I do wrong, what went wrong, why?
> with friendly greetings,
> P_ter
There is not much hope to working with the infinite summation since
Sum cannot compute a closed form for the result. Hence anything it
does will rely on some sort of approximation via numeric summing of
some initial summands, and extrapolation for the tail. This might or
might not give reliable results. Since it is fairly easy to bound the
error of finite summations "by hand", given that {p,q} are in the
whereabouts of {5,.01}, it is probably more reliable just to truncate
the sum at some sufficiently high, but finite, bound. 10000 will do
just fine (check it).
n[p_Real, q_Real] := Sum[k/(Exp[p + k*q] - 1), {k, 1, 10000}]
m[p_Real, q_Real] := Sum[1/(Exp[p + k*q] - 1), {k, 1, 10000}]
mm = 0.501;
mn = 41.0959;
In[13]:= rt =
FindRoot[{m[p, q] - mm, n[p, q] - mn}, {{p, 5}, {q, 0.1}}]
Out[13]= {p -> 5.09055, q -> 0.0122471}
In[14]:= m[p, q] /. rt
Out[14]= 0.501
In[15]:= n[p, q] /. rt
Out[15]= 41.0959
I'll sum more terms for n[] just to show that the result really does
not change (that is, the truncation approximation was fine).
In[26]:= Sum[1/(Exp[p + k q] - 1) /. rt, {k, 1, 10^5}] // InputForm
Out[26]//InputForm=
0.501000000000001
This is the same as I get if opnly summing to 10^4.
Daniel Lichtblau
Wolfram Research