Re: Re: How to get an answer as a Root object?
- To: mathgroup at smc.vnet.net
- Subject: [mg57189] Re: [mg57156] Re: [mg57137] How to get an answer as a Root object?
- From: DrBob <drbob at bigfoot.com>
- Date: Fri, 20 May 2005 04:43:23 -0400 (EDT)
- References: <200505170520.BAA25934@smc.vnet.net> <200505190708.DAA13114@smc.vnet.net>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
Two questions, Daniel:
1) Is this the missing definition of sroot?
sroot = FindRoot[expr == 0, {s, 2.7}]
2) Doesn't Drop do nothing at all in this statement?
fax2 = Drop[Map[First,FactorList[poly, Extension->Automatic]]];
Bobby
On Thu, 19 May 2005 03:08:33 -0400 (EDT), 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
>
>
>
>
>
--
DrBob at bigfoot.com
- 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?