MathGroup Archive 2004

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

Search the Archive

Re: Re: Re: Re: How to solve nonlinear equations?


DrBob wrote:
> Neither of us has answered the original question, which was "How to solve...".
> 
> Let's say we need 200 digits. How do you solve it?
> 
> Bobby
> [...]
>>"DrBob" <drbob at bigfoot.com> schrieb im Newsbeitrag
>>news:comgo1$89e$1 at smc.vnet.net...
>>
>>>The only trouble with THAT is that it doesn't work.
>>>
>>>eqns = {20.428836553499018 - Log[X1] == 468/67*X1 + 5790/1273*X2 -
>>>66257/1273*
>>>        X3 - 24660/1273*X4 - 79150/1273*X5, 17.011108423692498 -
>>>           Log[X2] == 5790/1273*X1 + 6294/1273*X2 - 66257/1273*X3 - 24660/
>>>          1273*X4 - 403.8069285*
>>>            X5, -29.72639373695347 - Log[X3] == -66257/1273*
>>>          X1 - 66257/1273*X2 - 2*
>>>        X3, -26.273726271581616 - Log[X4] == -24660/1273*X1 - 24660/
>>>          1273*X2 + 10.15330715*X4, -38.76695085346396 - Log[X5]
>>>== -79150/
>>>          1273*X1 - 403.8069285*X2 - 10.67374705*X5};
>>>
>>>FindRoot[eqns, {{X1, 1}, {X2, 1}, {X3, 1}, {X4, 1}, {X5, 1}}]
>>>eqns /. Equal -> Subtract /. %
>>>
>>>(FindRoot::lstol error)
>>>{X1 -> 0.5287215008048355 + 0.0078093799604902845*I,
>>>  X2 -> 0.002514442570050682 - 0.007506800892680369*I,
>>>  X3 -> 0.2750175495952261 + 0.0051666476248725625*I,
>>>  X4 -> -1.6157105189511614 -4.778719413301337*^-15*I,
>>>  X5 -> -0.006058532522555838 - 0.003759309847863656*I}
>>>{-0.0000172993518866571 - 5.25106242902848*^-7*I,
>>>  0.0014397436481736747 + 0.000060457451791506855*I,
>>>  -0.23588972784670026 + 0.007297504899863284*I,
>>>  -0.057825246363192695 + 3.1474540831359445*I,
>>>  0.0009303265262765592 + 0.00036856201666024546*I}
>>>
>>>Bobby
>>> [...]
>>>>
>>>>"Wei Wang" <weiwang at baosteel.com> schrieb im Newsbeitrag
>>>>news:cohj9d$1nr$1 at smc.vnet.net...
>>>>
>>>>>How to solve the following equations, where X1, X2, X3, X4 and X5 are
>>>>>variables?
>>>>>
>>>>>eqns = {20.428836553499018-ln(X1) ==
>>>>>468/67*X1+5790/1273*X2-66257/1273*X3-24660/1273*X4-79150/1273*X5,
>>>>>17.011108423692498-ln(X2) ==
>>>>>5790/1273*X1+6294/1273*X2-66257/1273*X3-24660/1273*X4-403.8069285*X5, -29.72639373695347-ln(X3)
>>>>>== -66257/1273*X1-66257/1273*X2-2*X3, -26.273726271581616-ln(X4)
>>>>>== -24660/1273*X1-24660/1273*X2+10.15330715*X4, -38.76695085346396-ln(X5)
>>>>>== -79150/1273*X1-403.8069285*X2-10.67374705*X5};
>>>>>
>>>
>>>--
>>>DrBob at bigfoot.com
>>>www.eclecticdreams.net

I'm not sure about the case of complex roots, but there is no real 
solution to this system. Showing that took me several steps though there 
may be a more enlightened method than what I reproduce below.

Start with

eqns = {20.428836553499018 - Log[X[1]] -
     (468/67*X[1] + 5790/1273*X[2] - 66257/1273*X[3] -
       24660/1273*X[4] - 79150/1273*X[5]),
   17.011108423692498 - Log[X[2]] -
     (5790/1273*X[1] + 6294/1273*X[2] - 66257/1273*X[3] -
	  24660/1273*X[4] - 403.8069285X[5]),
   -29.72639373695347 - Log[X[3]] -
     (-66257/1273*X[1] - 66257/1273*X[2] - 2*X[3]),
   -26.273726271581616 - Log[X[4]] -
     (-24660/1273*X[1] - 24660/1273*X[2] + 10.15330715*X[4]),
   -38.76695085346396 - Log[X[5]] -
     (-79150/1273*X[1] - 403.8069285*X[2] - 10.67374705*X[5])};

exprs = eqns /. Equal->Subtract;

First notice that we can replace each X[j] by its exponential. Moreover 
as all right hand sides are free of the logs, we know the X[j]'s must be 
positive. Hence we can use PowerExpand to in effect make those logs our 
new variables. I will also Rationalize though that's not strictly necessary.

newexprs = Rationalize[PowerExpand[exprs /. X[j_]:>Exp[lX[j]]], 0]

I now use some elimination ideals to obtain various expressions that 
must vanish. From these certain inequality deductions may be made.

expr34 = GroebnerBasis[newexprs, Exp[Array[lX,5]],
   Complement[Array[lX,5],{lX[3],lX[4]}]];

We first investigate what happens if lX[4] is nonnegative.

In[280]:= InputForm[e3 = First[expr34] /. lX[4]->Interval[{0,Infinity}]]
Out[280]//InputForm=
647238762444175200000000*E^lX[3] +
  Interval[{22053533331180859235803693, Infinity}] -
  323619381222087600000000*lX[3]

The next steps perhaps obvious to many people. We first see that this is 
strictly positive if lX[3] is negative.

In[289]:= e3 /. lX[3]->Interval[{-Infinity,0}]
Out[289]= Interval[{22053533331180859235803693, Infinity}]

We now check that the derivative with respect to lX[3] is strictly 
positive at for positive lX[3]. This guarantees that e3 cannot get 
smaller, hence is positive for all real lX[3].

In[290]:= D[e3,lX[3]] /. lX[3]->Interval[{0,Infinity}]
Out[290]= Interval[{323619381222087600000000, Infinity}]

We conclude at this point that we must require lX[4] to be negative. 
Taking this into account we will show that lX[1]<1 by showing that if we 
assume otherwise then a certain expression in our list will be strictly 
positive.

In[291]:= N[newexprs[[-2]] /. {lX[1]->Interval[{1,Infinity}],
   lX[4]->Interval[{-Infinity,0}],
   lX[2]->Interval[{-Infinity,Infinity}]}]
Out[291]= Interval[{16.2303, Infinity}]

Similarly we see we must have lX[2]<1.

In[293]:= N[newexprs[[-2]] /. {lX[2]->Interval[{1,Infinity}],
   lX[4]->Interval[{-Infinity,0}],
   lX[1]->Interval[{-Infinity,Infinity}]}]
Out[293]= Interval[{16.2303, Infinity}]

We now show these each must in fact be negative.

expr5 = First[GroebnerBasis[newexprs, Exp[Array[lX,5]],
   Complement[Array[lX,5],{lX[5]}]]];

e15 = expr5 /. {lX[1]->Interval[{0,1}], lX[2]->Interval[{-Infinity,1}]};

First show this is strictly positive for all negative lX[5].

In[297]:= e15 /. lX[5]->Interval[{-Infinity,0}]
Out[297]= Interval[{1022379647864000000, Infinity}]

Next show the derivative with respect to lX[5] is strictly positive for 
lX[5]>0.

In[299]:= D[e15,lX[5]] /. lX[5]->Interval[{0,Infinity}]
Out[299]= Interval[{422497213422049953, Infinity}]

This shows e15 cannot vanish for any real lX[5]. We can similarly show that

e25 = {lX[2]->Interval[{0,1}], lX[1]->Interval[{-Infinity,1}]};

cannot vanish for any lX[5]. We now know that {lX[1],lX[2],lX[4]} are 
all negative. We move to another elimination ideal.

expr3 = First[GroebnerBasis[newexprs, Exp[Array[lX,5]],
   Complement[Array[lX,5],{lX[3]}]]]

We see that this is strictly positive if lX[3] is 3.

In[311]:= (expr3 /. {lX[3]->3,
   lX[1]->Interval[{-Infinity,0}], lX[2]->Interval[{-Infinity,0}]}) > 0
Out[311]= True

The derivative with respect to lX[3] is always positive for lX[3]>3.

In[313]:= (D[expr3,lX[3]] /. lX[3]->Interval[{3,Infinity}]) > 0
Out[313]= True

Thus expr3 cannot vanish for lX[3]>=3, so lX[3]<3. We next look at one 
of our original expressions when lX[3]<-30. We will see it must in this 
case be strictly positive.

In[316]:= (newexprs[[3]] /. {lX[3]->Interval[{-Infinity,-30}],
   lX[2]->Interval[{-Infinity,0}], lX[1]->Interval[{-Infinity,0}]}) > 0
Out[316]= True

We now know -30<lX[3]<3. Finally we use all this in the first of newexprs.

In[318]:= (newexprs[[1]] /. {lX[3]->Interval[{-30,3}],
   lX[2]->Interval[{-Infinity,0}],
   lX[4]->Interval[{-Infinity,0}], lX[5]->Interval[{-Infinity,Infinity}],
   lX[1]->Interval[{-Infinity,0}]}) > 0
Out[318]= True

As it is strictly positive there are no real values that make newexprs 
simultaneously vanish.


Daniel Lichtblau
Wolfram Research



  • Prev by Date: Re: Re: A problem of numerical precision
  • Next by Date: How to evaluate Exp[I Pi(1+x)]?
  • Previous by thread: Re: Re: Re: How to solve nonlinear equations?
  • Next by thread: GUIKit and VisibleRect