Re: Re: How to get an answer as a Root object?
- To: mathgroup at smc.vnet.net
- Subject: [mg57195] Re: [mg57156] Re: [mg57137] How to get an answer as a Root object?
- From: Chris Chiasson <chris.chiasson at gmail.com>
- Date: Fri, 20 May 2005 04:43:36 -0400 (EDT)
- References: <200505170520.BAA25934@smc.vnet.net> <200505190708.DAA13114@smc.vnet.net>
- Reply-to: Chris Chiasson <chris.chiasson at gmail.com>
- Sender: owner-wri-mathgroup at wolfram.com
Daniel,
Thank you for taking the time to illustrate how to solve for
complicated roots using those exact methods.
In Method 2, how do you know which part of gb is the correct one?
In Method 1, there were some undefined variables, namely act and
sroot. I fixed the example code below:
srootval=s/.(sroot=
MapAt[Re,
FindRoot[expr,{s,2.7},WorkingPrecision->240,
PrecisionGoal->120],{1,2}]);
<<NumberTheory`Recognize`
fax=Drop[Map[First,FactorList[Recognize[srootval,30,s]]],1]
minpoly=First[Select[fax,Chop[#/.sroot]===0&]]
sexactroot=Solve[minpoly==0][[2]]
sexactroot//N
Thank you for your time,
On 5/19/05, Daniel Lichtblau <danl at wolfram.com> wrote:
> 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
>
>
>
--
Chris Chiasson
http://chrischiasson.com/
1 (810) 265-3161
- References:
- How to get an answer as a Root object?
- From: "David W. Cantrell" <DWCantrell@sigmaxi.org>
- Re: How to get an answer as a Root object?
- From: Daniel Lichtblau <danl@wolfram.com>
- How to get an answer as a Root object?