WOLFRAM

Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 1998

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Re: FindRoot accuracy/precision

Daniel Lichtblau (danl@wolfram.com) writes:
> Michael Chang wrote:
>> 
>> Thanks again for your response!
>> [--said response deleted for brevity-- dl]
>> 
>> Okay, I've had Solve running for two days now trying to solve my actual
>> seven variable problem, but haven't gotten any output yet.  (I'm
>> running a NEC P220MHz PC with Windoze 95 and 64MB RAM and Mathematica
>> 3.0.)
>> 
>> Specifically, my seven equations are as follows:
>> 
>> eqs={369664/981 - 1/(x2*(x3+x4))==0,
>> -6513786578938559/9223372036854775808+x1*x4+x1*x5+x6*x7==0,
>> -3858337547701687/9444732965739290427392+x1*x2*x4*x5+x1*x4*x6*x7+
>>         x1*x5*x6*x7==0,
>> -1/97280000000+x1*x2*x4*x5*x6*x7==0, 981/369664-x2*x3-x2*x4==0,
>> 77/924160000-x1*x2*x3*x4-x1*x2*x3*x5-x2*x3*x4*x6-x2*x3*x6*x7-
>>         x2*x4*x6*x7==0,
>> 1/1848320000000-x1*x2*x3*x4*x6*x7-x1*x2*x3*x5*x6*x7==0}
>> 
>> (The fractions above may seem ridiculous, but they give me infinite
>> precision.) So ... I try
>> 
>> ans=Solve[eqs,{x1,x2,x3,x4,x5,x6,x7}]
>> 
>> Thus far, I haven't had any output (unfortunately) after about 1.5 days.
>> Now I'm a newbie at Mathematica, so is there something that I'm doing
>> wrong here?  I initially started just wanted to find positive real root
>> solutions to these seven equations, but now I'm just struggling to find
>> *any* complex solution to this set of equations.  Sigh ...
>> 
>> Again, I'd be grateful for any suggestions or help that anyone may have!
>> 
>> Thanks again!
>> 
>> Mike
> 
> 
> Okay, I have some more suggestions. First, abort your Solve. It won't
> finish.
> 
> Here is what I tried, roughly in the order attempted. 
> 
> max = 50;
> start := Table[Random[Real,max,100], {Length[vars]}] Timing[rt =
> FindRoot[Evaluate[Thread[polys==0]],
> 	Evaluate[Sequence @@ Thread[{vars,start}]],
> 	WorkingPrecision->100, AccuracyGoal->20,
> 	MaxIterations->300]]
> 
> I got a failure-to-converge warning and a bad residual, so I upped the
> MaxIterations to 1000 (same story) and then 10000 (reprise). So I start
> to wonder about the existence of real solutions. Turns out they exist
> but are tricky to find this way. And for a good reason, to be
> disclosed.
> 
> Next I make explicit polynomials, in preparation for using
> GroebnerBasis.
> 
> p2 = Numerator[Together[polys]]
> 
> Now I compute a Groebner basis that in effect triangulates the system of
> polynomials. I first do this with a nonzero modulus, so as to keep
> coefficient growth from murdering the computation (it turns out this is
> not a problem, but often it is). We will see that we do not have seven
> independent polynomials, only six. What this means is any Newton's
> method based algorithm will have trouble because we do not have a
> nonsingular Jacobian; this explains why FindRoot was never successful.
> 
> In[29]:= Timing[gb = GroebnerBasis[p2, vars, Modulus->Prime[1111]];]
> Out[29]= {0.27 Second, Null}
> 
> In[30]:= InputForm[gb]
> Out[30]//InputForm=
> {1966 + 1942*x6*x7 + 5511*x6^2*x7^2 + 5810*x6^3*x7^3,
>  6312*x5 + 1894*x7 + 7949*x6*x7^2 + 339*x6^2*x7^3,
>  3020*x4 + 5130*x7 + 923*x6*x7^2 + 2114*x6^2*x7^3,
>  2507*x3 + 72*x7 + 3082*x6*x7^2 + 566*x6^2*x7^3,
>  809*x2 + 2999*x6 + 4901*x6^2*x7 + 5360*x6^3*x7^2,
>  957*x1 + 1371*x6 + 4793*x6^2*x7 + 6656*x6^3*x7^2}
> 
> Note that the first polynomial is not in x7 alone, but has x6 as well.
> Indeed, fortuitiously, it can be written as a polynomial in x6*x7.
> Before preceding further we will recompute this in characteristic zero.
> 
> In[31]:= Timing[gb = GroebnerBasis[p2, vars];] Out[31]= {0.72 Second,
> Null}
> 
> Now make that first polynomial univariate by algebraically replacing
> x6*x7 by z.
> 
> In[33]:=  InputForm[pol1 =
> 	Last[PolynomialReduce[gb[[1]], x6*x7 - z, {x6,x7,z}]];]
> Out[33]//InputForm=
> -144115188075855872 + 5727219797369691640625*z -
>  9900955599986609680000000*z^2 + 14019525496019259228160000000*z^3
> 
> Now let us find the roots.
> 
> In[35]:= InputForm[rz = NRoots[pol1==0, z, 50]] Out[35]//InputForm=
> z == 0.00002631578947368421257295093734440182212051317541593`49.3587 ||
>  z == 0.0003399551832627632631436772994454357857827081508994`49.4544 -
>    0.0005244573132036400673608051786757841450315497819155`49.6427*I ||
>  z == 0.0003399551832627632631436772994454357857827081508994`49.4544 +
>    0.0005244573132036400673608051786757841450315497819155`49.6427*I
> 
> We have a positive root, which means we might have solutions for the
> original system for which x6 and x7 are positive, and their product
> must then be approximately 0.0000263158. We can now play with FindRoot
> a bit more, giving these new polynomials and adjoining a specific value
> for x7, say, to have #equations = #unknowns.
> 
> By trial and error I found I had bad luck with MaxIterations set to 100,
> but better results at 300.
> 
> In[49]:= Timing[rt = FindRoot[Evaluate[Thread[Prepend[gb,x7-1]==0]],
>         Evaluate[Sequence @@ Thread[{vars,start}]],
>         WorkingPrecision->100, AccuracyGoal->20,
>         MaxIterations->300];]
> 
> Out[49]= {2.58 Second, Null}
> In[50]:= polys /. rt
>                -97    PolynomialReduce[gb, x6*x7 -  0.0000263158]]  
> -103       -106       -110       -102       -106
> Out[50]= {0. 10   , 0. 10    , 0. 10    , 0. 10    , 0. 10    , 0. 10   
> ,
>           -111
>>    0. 10    }
> 
> Real roots, and small residuals to boot. But the roots are not all
> positive. By looking hard at gb, and keeping in mind what the product
> x6*x7 must be, we can prove that there is no nonnegative solution to
> the system.
> 
> In[62]:= Map[Last, PolynomialReduce[gb, x6*x7 - rz[[1,2]]]] // InputForm
> Out[62]//InputForm=
> {3.239236523077941139396279002864633272546781779982`49.3587*^-33,
>  1.0672319074126148445599036240762780763974738580124`49.3587*^86*x5 +
>   2.661030840495747972280200538958840019865272745528`49.3587*^86*x7,
>  1.329227995784915872903807060280344576`49.3587*^36*x4 +
>   1.1820175769042802462402839996427728924163754699203`49.3587*^37*x7,
>  3.483354311655665952710048388609092358440802231132`49.3587*^87*x3 +
>   3.570195340839940700543429000179020127358523417341`49.3587*^86*x7,
>  7.714912756920962485048802864841746175066733934288`49.3587*^124*x2 +
>   8.649182463653030159385948370295392422188602313057`49.3587*^125*x6,
>  2.466352649960385217966927145432296886685259963593`49.3587*^100*x1 +
>   5.596581083779066501874325814280057485830651652433`49.3587*^100*x6}
> 
> As expected, we get a small value for the first polynomial, because it
> is supposed to vanish for that value of x6*x7. The rest give
> polynomials that are each linear in x7 and in one other variable, in
> effect solving for that variable in terms of x7. We see, moreover, that
> all coefficients are positive. So after dividing each polymonial by its
> x7 coefficient we have eqns of the form
> 
> positive*xj+x7==0.
> 
> This gives a negative solution for each xj for 1<=j<=5. Conclusion as
> stated above: there is no root of the original system with
> positive-only values. Which agrees with my experience using FindRoot on
> gb several times in the manner shown above.
> 
> 
> Daniel Lichtblau
> Wolfram Research
>
  • Prev by Date: Mathematica 3.0: How to change thickness of 3D-curve (ParametricPlot3D)
  • Next by Date: RE: commutativity
  • Prev by thread: Re: Re: Re: FindRoot accuracy/precision
  • Next by thread: Re: FindRoot accuracy/precision