# Re: Re: Re: FindRoot accuracy/precision

*To*: mathgroup@smc.vnet.net*Subject*: [mg11506] Re: [mg11460] Re: Re: FindRoot accuracy/precision*From*: Daniel Lichtblau <danl@wolfram.com>*Date*: Fri, 13 Mar 1998 12:21:58 -0500*References*: <199803030411.XAA02345@smc.vnet.net.> <6dk5kd$oam$3@dragonfly.wolfram.com> <199803060541.AAA22147@smc.vnet.net.> <6dtb42$s4u$3@dragonfly.wolfram.com> <199803120634.BAA23619@smc.vnet.net.>

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

**References**:**FindRoot accuracy/precision***From:*bt585@FreeNet.Carleton.CA (Michael Chang)

**Re: FindRoot accuracy/precision***From:*bt585@FreeNet.Carleton.CA (Michael Chang)

**Re: Re: FindRoot accuracy/precision***From:*bt585@FreeNet.Carleton.CA (Michael Chang)

**Re: MultipleListPlot with Log Linear Scale**

**Re: Question about Coefficient**

**Re: Re: FindRoot accuracy/precision**

**Re: Re: Re: FindRoot accuracy/precision**