Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Cannot NSolve a system of equations

• To: mathgroup at smc.vnet.net
• Subject: [mg88945] Re: Cannot NSolve a system of equations
• From: Szabolcs Horvát <szhorvat at gmail.com>
• Date: Thu, 22 May 2008 02:34:51 -0400 (EDT)
• Organization: University of Bergen
• References: <200805201052.GAA05057@smc.vnet.net> <g11r0h$a7v$1@smc.vnet.net>

Daniel Lichtblau wrote:
> murat.koyuncu at gmail.com wrote:
>> Dear all,
>>
>> I have the following system that I need to solve, but I cannot get a
>> sensible result.
>>
>> Unprotect[In,Out];Clear[In,Out];ClearAll["Global`*"];
>> zet=0.083;
>> phi=0.75;eta=1.75;alpha=0.64;y1=0.235457064;y2=0.512465374;y3=0.781779009;
>> y4=1.109572176; y5=2.360726377;tau1=zet y1^phi;tau2=zet
>> y2^phi;tau3=zet y3^phi;tau4=zet y4^phi;tau5=zet y5^phi;
>>  taubar=(tau1 y1+tau2 y2+tau3 y3+tau4 y4+tau5 y5)/
>> 5;a1=(1+phi)tau1;a2=(1+phi)tau2;a3=(1+phi)tau3;a4=(1+phi)tau4;a5=(1+phi)tau5;
>>
>> eqns1={x1==(roverw(1-tau1)+((roverw+(1-x))y1-1)( (1-taubar+(1-abar)/
>> eta)x+(taubar-tau1)roverw-(1-taubar)))/(roverw (1-tau1+(1-a1)/eta)-
>> ( (1-taubar+(1-abar)/eta)x+(taubar-tau1)roverw-(1-taubar))),
>> x2==(roverw(1-tau2)+((roverw+(1-x))y2-1)( (1-taubar+(1-abar)/eta)x+
>> (taubar-tau2)roverw-(1-taubar)))/(roverw (1-tau2+(1-a2)/eta)-( (1-
>> taubar+(1-abar)/eta)x+(taubar-tau2)roverw-(1-taubar))),
>> x3==(roverw(1-tau3)+((roverw+(1-x))y3-1)( (1-taubar+(1-abar)/eta)x+
>> (taubar-tau3)roverw-(1-taubar)))/(roverw (1-tau3+(1-a3)/eta)-( (1-
>> taubar+(1-abar)/eta)x+(taubar-tau3)roverw-(1-taubar))),
>> x4==(roverw(1-tau4)+((roverw+(1-x))y4-1)( (1-taubar+(1-abar)/eta)x+
>> (taubar-tau4)roverw-(1-taubar)))/(roverw (1-tau4+(1-a4)/eta)-( (1-
>> taubar+(1-abar)/eta)x+(taubar-tau4)roverw-(1-taubar))),
>> x5==(roverw(1-tau5)+((roverw+(1-x))y5-1)( (1-taubar+(1-abar)/eta)x+
>> (taubar-tau5)roverw-(1-taubar)))/(roverw (1-tau5+(1-a5)/eta)-( (1-
>> taubar+(1-abar)/eta)x+(taubar-tau5)roverw-(1-taubar))),
>> x==(x1+x2+x3+x4+x5)/5, abar == (a1 x1+a2 x2+a3 x3+a4 x4+a5 x5)/
>> (5x),roverw==(1-x)(1-alpha)/alpha };
>>
>> sol=NSolve[eqns1,{x, x1,x2,x3,x4,x5,abar, roverw}];
>>
>> eqns1 /. sol
>>
>> Out[741]={{False, False, False, False, False, True, True, True},
>> {False, False,
>>    False, False, False, True, True, True}, {False, False, False,
>>   False, False, True, False, True}, {False, False, True, False, False,
>>    True, False, True}, {False, False, True, False, False, True, False,
>>    True}, {False, False, True, True, False, True, True, False}}
>>
>>
>> What am I doing wrong? Is it just because the system is too
>> complicated?
>>
>> Any help would be truly appreciated.
>> Murat
> 
> Quite possibly there are issues involving numeric stability and the 
> presence of denominators. I was able to get a solution set, containg two 
> solutions, by starting with exact input and then numericizing to high 
> precision.
> 
> zet = 83/1000;
> phi = 3/4;
> eta = 7/4;
> alpha = 16/25;
> 
> {y1,y2,y3,y4,y5} = Rationalize[
>    {0.235457064,0.512465374,0.781779009,1.109572176,2.360726377}, 0];
> 
> With this I can do:
> 
> In[34]:=  Timing[sol = NSolve[N[eqns,500],vars];]
> Out[34]= {5., Null}
> 
> In[36]:= eqns/.sol
> Out[36]= {{True, True, True, True, True, True, True, True},
>     {True, True, True, True, True, True, True, True}}
> 
> Here are the solution values, at machine precision.
> 
> In[39]:= InputForm[N[sol]]
> Out[39]//InputForm=
> {{x -> 0.6251836373550925, x1 -> 0.6454246796060954,
>    x2 -> 0.6474758190687543, x3 -> 0.6445628507802813,
>    x4 -> 0.6355486246954559, x5 -> 0.5529062126248758,
>    abar -> 0.13411739130856176, roverw -> 0.21083420398776043},
>   {x -> 1., x1 -> 1., x2 -> 1., x3 -> 1., x4 -> 1., x5 -> 1.,
>    abar -> 0.13829898904658144, roverw -> 0.}}
> 

An equivalent, but somewhat simpler approach is to simply increase the 
nominal precision of the numbers:

In[25]:= eqn2 = SetPrecision[eqns1, 30];

In[26]:= sol = NSolve[eqn2, {x, x1, x2, x3, x4, x5, abar, roverw}];

In[27]:= eqn2 /. Equal -> Subtract /. sol

Out[27]= {{0.*10^-25, 0.*10^-25, 0.*10^-23, 0.*10^-24, 0.*10^-25,
   0.*10^-27, 0.*10^-28, 0.*10^-28}, {0.*10^-24, 0.*10^-24, 0.*10^-24,
   0.*10^-23, 0.*10^-23, 0.*10^-26, 0.*10^-27,
   0.*10^-27}, {0.*10^-25 + 0.*10^-25 I, 0.*10^-26 + 0.*10^-26 I,
   0.*10^-26 + 0.*10^-26 I, 0.*10^-26 + 0.*10^-26 I,
   0.*10^-26 + 0.*10^-26 I, 0.*10^-28 + 0.*10^-28 I,
   0.*10^-28 + 0.*10^-28 I,
   0.*10^-29 + 0.*10^-29 I}, {0.*10^-25 + 0.*10^-25 I,
   0.*10^-26 + 0.*10^-26 I, 0.*10^-26 + 0.*10^-26 I,
   0.*10^-26 + 0.*10^-26 I, 0.*10^-26 + 0.*10^-26 I,
   0.*10^-28 + 0.*10^-28 I, 0.*10^-28 + 0.*10^-28 I,
   0.*10^-29 + 0.*10^-29 I}, {0.*10^-29, 0.*10^-29, 0.*10^-29,
   0.*10^-29, 0.*10^-29, 0.*10^-29, 0.*10^-30, 0.*10^-29}, {0.*10^-28,
   0.*10^-28, 0.*10^-28, 0.*10^-28, 0.*10^-28, 0.*10^-29, 0.*10^-29,
   0.*10^-29}}

• References:
  • Cannot NSolve a system of equations
    • From: murat.koyuncu@gmail.com