Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Cannot NSolve a system of equations

• To: mathgroup at smc.vnet.net
• Subject: [mg88906] Re: Cannot NSolve a system of equations
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 21 May 2008 14:48:20 -0400 (EDT)
• References: <200805201052.GAA05057@smc.vnet.net>

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.}}


Daniel Lichtblau
Wolfram Research

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