[Date Index]
[Thread Index]
[Author Index]
Re: How to get an answer as a Root object?
*To*: mathgroup at smc.vnet.net
*Subject*: [mg57156] Re: [mg57137] How to get an answer as a Root object?
*From*: Daniel Lichtblau <danl at wolfram.com>
*Date*: Thu, 19 May 2005 03:08:33 -0400 (EDT)
*References*: <200505170520.BAA25934@smc.vnet.net>
*Sender*: owner-wri-mathgroup at wolfram.com
David W. Cantrell wrote:
> This question must have an obvious answer!
>
> A certain expression, given below my signature, involves nothing worse than
> some nested square roots. It equals 0 for a value of the variable s near
> 2.707, easily approximated using FindRoot. However, I would like to get
> that value expressed precisely as a Root object. Trying to use Solve or
> Reduce to do that, Mathematica seems to take "forever".
>
> Of course, I could do this by hand; in the days before CASs -- before I got
> spoiled by them! -- I wouldn't have hesitated to do so. But surely there
> must be an easy way to get Mathematica to give that value as a Root object.
> How?
>
> David Cantrell
>
>
> (-1 + Sqrt[2] - Sqrt[(1 + 2*Sqrt[2] - s)*(-1 + s)] + s)*(-1 - Sqrt[2] +
> Sqrt[2]*s - Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] - s)*s]) -
> (-1 - Sqrt[2] + Sqrt[2]*s + Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] - s)*s])*
> (-2 - 2*Sqrt[2] + 2*s + 2*Sqrt[1 - (1/4)*(-1 - Sqrt[2] + Sqrt[(1 +
> 2*Sqrt[2] - s)*(-1 + s)] + s)^2])
There are a few reasonable methods. Here are some.
expr = (-1 + Sqrt[2] - Sqrt[(1 + 2*Sqrt[2] - s)*(-1 + s)] + s)*(-1 -
Sqrt[2] +
Sqrt[2]*s - Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] - s)*s]) -
(-1 - Sqrt[2] + Sqrt[2]*s + Sqrt[-1 - 2*Sqrt[2] + 2*(2 + Sqrt[2] -
s)*s])*
(-2 - 2*Sqrt[2] + 2*s + 2*Sqrt[1 - (1/4)*(-1 - Sqrt[2] + Sqrt[(1 +
2*Sqrt[2] - s)*(-1 + s)] + s)^2]);
Method 1.
First we use numerical approximation to recognize a root polynomial.
Due, I suspect, to nested radicals, we'll get a small spurious imaginary
part which I remove using Re.
In[5]:= InputForm[srootval = Re[s/.sroot]]
Out[5]//InputForm=
2.706945837045705738269073202812169790793196895164013983001575249106155514352\
266542032260316431664549796229819466076004645196713561849037`120.
To find a candidate polynomial we'll use Recognize (uses a method based
on lattice reduction), then factor, and find the factor with the right root.
<<NumberTheory`Recognize`
fax = Drop[Map[First,FactorList[Recognize[srootval, 30, s]]],1];
InputForm[minpoly = First[Select[fax, Chop[#/.sroot,10^(-acc)]===0&]]]
Why is the factorization needed? Because the larger product polynomial
has smaller vector norm in the lattice than does the actual minimal
polynomial.
Method 2.
We'll compute a Groebner basis with 's' the main variable. The hope,
born out in this case, is that internal handling of algebraics will
force to the fore a "polynomial" in s.
gb = GroebnerBasis[expr, s];
I do not show it here, but the polynomial we want is readily seen to be
the first one. Again we'll factor and find the one with the correct
root. As there are algebraics in the coefficients of this polynomial we
use Extension->Automatic in the factoring.
poly = First[gb];
fax2 = Drop[Map[First,FactorList[poly, Extension->Automatic]]];
In[21]:= InputForm[minpoly2 = First[Select[fax2,
Chop[#/.sroot,10^(-acc)]===0&]]]
Out[21]//InputForm=
-136 - 96*Sqrt[2] + (512 + 364*Sqrt[2])*s + (-804 - 578*Sqrt[2])*s^2 +
(676 + 486*Sqrt[2])*s^3 + (-329 - 220*Sqrt[2])*s^4 + (96 +
44*Sqrt[2])*s^5 -
16*s^6
How to represent as a polynomial over the integers? Just multiply by an
algebraic conjugate.
sqrtConjugate[ee_] := ee /. aa_^(1/2)->-aa^(1/2)
In[23]:= InputForm[minpoly2 = Expand[minpoly2*sqrtConjugate[minpoly2]]]
Out[23]//InputForm=
64 + 512*s - 6112*s^2 + 21024*s^3 - 32136*s^4 + 10832*s^5 + 43568*s^6 -
86152*s^7 + 81425*s^8 - 46080*s^9 + 15872*s^10 - 3072*s^11 + 256*s^12
Method 3.
Replace radicals by hand, use defining relations explicitly, formulate
Groebner basis, find polynomial in s alone.
rep1 = {Sqrt[2]->rad0};
rep2 = {Sqrt[(1+2*Sqrt[2]-s)*(-1+s)] -> rad1} /. rep1;
rep3 = ({Sqrt[-1-2*Sqrt[2]+2*(2+Sqrt[2]-s)*s] -> rad2,
Sqrt[1-(1/4)*(-1-Sqrt[2]+Sqrt[(1+2*Sqrt[2]-s)*(-1+s)]+s)^2] ->
rad3} /. rep1 /. rep2)
repssq = Map[#[[1]]^2-#[[2]]^2&, Flatten[{rep1,rep2,rep3}]]
expr2 = (expr /. rep1) /. reps;
In[64]:= InputForm[Timing[poly4 =
First[GroebnerBasis[Join[repssq,{expr2}], s,
{rad0,rad1,rad2,rad3}, MonomialOrder->EliminationOrder]]]]
Out[64]//InputForm=
{0.30000000000000115*Second, 200704 + 3956736*s - 16568320*s^2 -
259776512*s^3 + 2455161856*s^4 - 9217882112*s^5 + 16736476672*s^6 -
4378572288*s^7 - 50113953728*s^8 + 119132812800*s^9 - 109419876224*s^10 -
38011997376*s^11 + 225422464224*s^12 - 259623250704*s^13 +
80225006088*s^14 + 142924054128*s^15 - 196411664327*s^16 +
53713807776*s^17 + 132853129136*s^18 - 220218741248*s^19 +
191796620928*s^20 - 115209149440*s^21 + 51080841728*s^22 -
17021165568*s^23 + 4236362752*s^24 - 766771200*s^25 + 95551488*s^26 -
7340032*s^27 + 262144*s^28}
Again one would factor and find the one with the correct root.
I'm sure there is a resultant formulation as well that will get a
polynomial in s for the root in question.
Daniel Lichtblau
Wolfram Research
Prev by Date:
**Re: UpSetDelayed and N**
Next by Date:
**Errors from FindFit**
Previous by thread:
**How to get an answer as a Root object?**
Next by thread:
**Re: Re: How to get an answer as a Root object?**
| |