Prev: ReplaceAll behavour
Next: Mathematica 6 and system of matrix and scalar ODEs in NDSolve:
From: P_ter on 30 Oct 2007 04:40
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
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

From: Jean-Marc Gulliet on 31 Oct 2007 07:16
P_ter wrote:

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

Mixing symbolic and numeric result might not be the best approach. If,
say, Sum[1/(Exp[p + k q] - 1), {k, 1, Infinity}] has no closed form (as
Mathematica claims), you are doomed to fail with this approach.

You may be more successful by using purely numerical methods and
controlling the precision as in the following example (though I do not
know whether the result {p -> 3.74867, q -> 0.0464572} makes sense).

In[1]:= Block[{$MinPrecision = 200, $MaxPrecision = 200},
n[p_?NumericQ, q_?NumericQ] :=
NSum[k/(Exp[p + k q] - 1), {k, 1, 10000}];
m[p_?NumericQ, q_?NumericQ] :=
NSum[1/(Exp[p + k q] - 1), {k, 1, 10000000}];
mm = 0.501;
mn = 41.0959;
FindRoot[{m[p, q] - mm, n[p, q] - mn}, {{p, 5}, {q, 0.1}}]
]

During evaluation of In[1]:= NIntegrate::slwcon: Numerical \
integration converging too slowly; suspect one of the following: \
singularity, value of the integration is 0, highly oscillatory \
integrand, or WorkingPrecision too small. >>

During evaluation of In[1]:= NIntegrate::ncvb: NIntegrate failed to \
converge to prescribed accuracy after 9 recursive bisections in k \
near {k} = {120.855}. NIntegrate obtained -4.98592*10^7 and \
136833.21646408358` for the integral and error estimates. >>

During evaluation of In[1]:= SequenceLimit::seqlim: The general form \
of the sequence could not be determined, and the result may be \
incorrect. >>

During evaluation of In[1]:= NIntegrate::slwcon: Numerical \
integration converging too slowly; suspect one of the following: \
singularity, value of the integration is 0, highly oscillatory \
integrand, or WorkingPrecision too small. >>

During evaluation of In[1]:= NIntegrate::ncvb: NIntegrate failed to \
converge to prescribed accuracy after 9 recursive bisections in k \
near {k} = {74.3418}. NIntegrate obtained -4.99961*10^7 and \
2600.0866205544694` for the integral and error estimates. >>

During evaluation of In[1]:= NIntegrate::slwcon: Numerical \
integration converging too slowly; suspect one of the following: \
singularity, value of the integration is 0, highly oscillatory \
integrand, or WorkingPrecision too small. >>

During evaluation of In[1]:= General::stop: Further output of \
NIntegrate::slwcon will be suppressed during this calculation. >>

During evaluation of In[1]:= NIntegrate::ncvb: NIntegrate failed to \
converge to prescribed accuracy after 9 recursive bisections in k \
near {k} = {231.1}. NIntegrate obtained -4.99785*10^7 and \
24827.157779411496` for the integral and error estimates. >>

During evaluation of In[1]:= General::stop: Further output of \
NIntegrate::ncvb will be suppressed during this calculation. >>

During evaluation of In[1]:= SequenceLimit::seqlim: The general form \
of the sequence could not be determined, and the result may be \
incorrect. >>

During evaluation of In[1]:= SequenceLimit::seqlim: The general form \
of the sequence could not be determined, and the result may be \
incorrect. >>

During evaluation of In[1]:= General::stop: Further output of \
SequenceLimit::seqlim will be suppressed during this calculation. >>

Out[1]= {p -> 3.74867, q -> 0.0464572}

Regards,
--
Jean-Marc

From: P_ter on 1 Nov 2007 06:51
Hello,
I checked the sums for stability. Seems ok up until i=100000.
Also: I did Newton-Raphson to solve this problem, into quite a few decimals.

But, I could not do it with Mathematica's FindRoot.
The two values mm and mn result in a first(!!) order approximation to p and q given in my mail. And reverse.

It is that I do not understand why such a beast as FindRoot comes with so nasty messages. Of course I accept that finding a root has a kind of art in it.
But calculating myself a first order approximation in two simple lines of code with only Log[] in it and that FindRoot can not find the roots, still urges me to ask: what did I do wrong?
with friendly greetings,
P_ter

From: danl on 2 Nov 2007 04:39
On Oct 30, 3:40 am, P_ter <peter_van_summe...(a)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




 | 
Pages: 1
Prev: ReplaceAll behavour
Next: Mathematica 6 and system of matrix and scalar ODEs in NDSolve: