Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 1998

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

Search the Archive

Re: Re: Re: FindRoot accuracy/precision

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: Re: MultipleListPlot with Log Linear Scale
  • Next by Date: Re: Question about Coefficient
  • Prev by thread: Re: Re: FindRoot accuracy/precision
  • Next by thread: Re: Re: Re: FindRoot accuracy/precision