Re: Re: How to get an answer as a Root object?
- To: mathgroup at smc.vnet.net
- Subject: [mg57198] Re: [mg57156] Re: [mg57137] How to get an answer as a Root object?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 20 May 2005 04:43:51 -0400 (EDT)
- References: <200505170520.BAA25934@smc.vnet.net> <200505190708.DAA13114@smc.vnet.net> <acbec1a405051916553ad63c12@mail.gmail.com>
- Sender: owner-wri-mathgroup at wolfram.com
Chris Chiasson wrote:
> 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,
See bottom. For better or worse I decided to keep the whole post so it
would at least be self contained.
> 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
I inadvertantly left out the actual FindRoot. It was something along the
lines you have. Here are the first few lines I managed to omit from my
last post.
acc = 100;
sroot = FindRoot[expr==0, {s,2.7}, AccuracyGoal->acc,
PrecisionGoal->acc, WorkingPrecision->(1.2*acc)];
As for figuring out which polynomial to use in the Groebner basis from
method 2, I knew because I looked. I should have stated this outright so
it did not resemble a magic trick. In brief, the expectation should be
that numeric algebraics e.g Sqrt[2] will be ordered lowest, so the
stated variable 's' will be interms of them. But algebraics containing s
will be ordered higher so that we do not get a circularity expressing s
in terms of an algebraic in s (which, in essence, is what the given
equation is doing).
In any case I would neither expect most users to realize this, not would
I necessarily trust it to work that way myself without checking. Which
is what I did.
To finish off by extracting a Root object, one can do as below. As ever,
we will select based on numerical value being equal to that we already
found, srootval.
In[21]:= InputForm[root = First[
Select[s/.Solve[minpoly2==0,s], #==srootval&]]]
Out[21]//InputForm=
1 + Sqrt[1 - Root[64 + 256*#1 - 1464*#1^2 + 1896*#1^3 -
495*#1^4 - 512*#1^5 + 256*#1^6 & , 1, 0]]
If we prefer a single root form, we can do
In[24]:= InputForm[root = RootReduce[root]]
Out[24]//InputForm=
Root[64 + 512*#1 - 6112*#1^2 + 21024*#1^3 - 32136*#1^4 + 10832*#1^5 +
43568*#1^6 - 86152*#1^7 + 81425*#1^8 - 46080*#1^9 + 15872*#1^10 -
3072*#1^11 + 256*#1^12 & , 8, 0]
My apologies for the omissions (which I hope this note will fix but I'm
not going to guarantee it).
Daniel Lichtblau
Wolfram Research
- Follow-Ups:
- Re: Re: Re: How to get an answer as a Root object?
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Re: Re: Re: How to get an answer as a Root object?
- 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?