Re: Solve for positive, real-valued solutions
- To: mathgroup@smc.vnet.net
- Subject: [mg12494] Re: [mg12415] Solve for positive, real-valued solutions
- From: "Jrgen Tischer" <jtischer@pitagoras.univalle.edu.co>
- Date: Tue, 19 May 1998 13:31:56 -0400
Hi Walter,
I'm concentrating on the long version.
First it's easy to see that there is for every c exactly one positive
solution (say by using Decartes rule). So to avoid having Mathematica
searching for all solutions one can search just for the first Root
which is positive. Now the Root objects are sorted, first coming the
real roots in natural order and afterwards the complex ones, so
r[c_, n_] := Root[#1^c/(c!*Sum[#1^i/i!, {i, 0, c}]) - 0.02 & , n]
posRoot[c_]:= Module[{sol,n=1}, While[(sol=r[c,n++])<=0];sol]
Now if you use this function, it works ok for not too big c, that is on
my PC it worked up to c==74 and then it came up with an awfully false
solution (something like 8 instead of 60). Of course this comes from
using a 0.02 together with real high powers. Now one can redefine the
above functions using 1/50 instead of .02 and then asking Mathematica
for more than 16 digits, so it will use higher precision and control
the error. This (of course) solves the problem of the false answers,
but is very slow. So to speed up the calculus I propose the following:
If one looks at the solutions for the first 50 or so c's, it becomes
obvious that the difference between two adjacent solutions is slowly
growing (looks like approaching 1 from below), so if one calculates
solutions, the last two solutions allow a rather good guess for the
next one. Now use FindRoot with this guess and an extra
WorkingPrecision and it will give fast and quite accurate answers ( I
calculated the first 250 in 24 seconds on a Pentium Pro 200).
In[1]:= solutions={r[1,1],r[2,2]}
Out[1]= {0.0204082,0.223467}
In[2]:=
nextSolution := Module[{guess = 2*Last[solutions] - solutions[[-2]],
c = Length[solutions] + 1,
sol},
sol = a /. FindRoot[50*a^c E^(-a) == Gamma[1
+ c, a], {a, guess},
WorkingPrecision -> 25];
AppendTo[solutions, sol];
sol]
In[3]:= Table[nextSolution,{250}];
This function in FindRoot was the result of using FullSimplify on your
function, I expect that makes it faster. (In Root you can't use it.)
Jürgen
-----Original Message-----
From: Walter Johnston <walterj@omega.uta.edu> To: mathgroup@smc.vnet.net
Subject: [mg12494] [mg12415] Solve for positive, real-valued solutions
>How can I tell Solve[] to give me only positive, real-valued solutions?
>
>Long Version:
>
>I am trying to build a table of values for use in a database. The
>expression is
>
>P[c_,a_]:=(a^c/c!)/Sum[a^i/i!,{i,0,c}].
>
>and I then tell the system to
>
>Do[Print[cc," => ",Solve[P[c,a]==0.02,a]],{cc,51,100,1}]
>
>which works, but gives me a lot more than I want.
>
>I have just started using Version 3.1 on my notebook (Windows NT 4.0)
>and would really appreciate any suggestions on how to speed this up and
>just get the positive, real-valued solutions I need for the database.
>
>(RTFM statements won't help much. I am reading, but trying to shorten
>the learning curve and get some solutions NOW.)
>
>Thanks.
>
>Walter Johnston
>
>