Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Cannot NSolve a system of equations

• To: mathgroup at smc.vnet.net
• Subject: [mg88937] Re: Cannot NSolve a system of equations
• From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
• Date: Wed, 21 May 2008 14:54:09 -0400 (EDT)
• References: <200805201052.GAA05057@smc.vnet.net>

On 20 May 2008, at 19:52, 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
>


You are doing nothing wrong, just demanding to much. Instead trying to  
verify the equations just calculate the residues:

eqns1 /. Equal -> Subtract /. sol
  {{-7.815970093361102*^-12, -8.151346264639869*^-11,
      -7.526068657170981*^-11, -2.701199264265597*^-10,
      -5.07498043589294*^-10, -9.947598300641403*^-14,
      8.881784197001252*^-16, 0.},
    {4.4158454670650826*^-11, 1.5007728393356956*^-9,
      3.164917217191032*^-10, -1.6515500078639889*^-10,
      -2.8396840434652404*^-11, -3.3306690738754696*^-16,
      -1.5265566588595902*^-15, 0.},
    {-4.1599612643494766*^-11 + 5.3216098194752703*^-11*
          I, -9.377387755193922*^-12 -
        1.566604623803869*^-10*I, 8.095479842040731*^-11 -
        1.2360334977756793*^-10*I,
      -7.734213269827706*^-11 + 7.310063665499911*^-11*I,
      2.5863755581667647*^-11 + 2.5732305175552028*^-11*
          I, 3.3306690738754696*^-16 -
        8.604228440844963*^-16*I,
      -2.525757381022231*^-15 + 8.770761894538737*^-15*I,
      2.7755575615628914*^-17 + 0.*I},
    {-4.1599612643494766*^-11 - 5.3216098194752703*^-11*
          I, -9.377387755193922*^-12 +
        1.566604623803869*^-10*I, 8.095479842040731*^-11 +
        1.2360334977756793*^-10*I,
      -7.734213269827706*^-11 - 7.310063665499911*^-11*I,
      2.5863755581667647*^-11 - 2.5732305175552028*^-11*
          I, 3.3306690738754696*^-16 +
        8.604228440844963*^-16*I,
      -2.525757381022231*^-15 - 8.770761894538737*^-15*I,
      2.7755575615628914*^-17 + 0.*I},
    {1.7763568394002505*^-15, -3.6415315207705135*^-14,
      4.052314039881821*^-15, 1.867950238931826*^-14,
      1.587618925213974*^-14, -5.551115123125783*^-17,
      1.2483070133129104*^-14, 0.},
    {1.5765166949677223*^-14, -8.182343691487404*^-14,
      1.2656542480726785*^-14, 3.6193270602780103*^-14,
      1.7430501486614958*^-14, 0.,
      -3.4139358007223564*^-15, -2.0358714714063808*^-14}}

You see that all the residues are small, so the fact that some of them  
are too large to be regarded by Mathematica as 0 is due to inexact  
coefficients and round off errors. For example:

Chop[%, 10^(-8)]
{{0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0},
    {0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0},
    {0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0}}

which means that the residues are all less than 10^-8, which seems  
pretty good to me. There is nothing more you can do with machine  
precision coefficients. A system with inexact coefficients does not,  
of course, have an "exact solution" in any meaninful sense. One thing  
you can do is rationalize your system and then compute an approximate  
solution to some degree of precision.


eqns2 = Rationalize[eqns1, 0];

sol2 = NSolve[eqns2, {x, x1, x2, x3, x4, x5, abar, roverw},  
WorkingPrecision -> 20];

Now of course:

eqns2 /. sol1
{{True, True, True, True, True, True, True, True},
    {True, True, True, True, True, True, True, True},
    {True, True, True, True, True, True, True, True},
    {True, True, True, True, True, True, True, True},
    {True, True, True, True, True, True, True, True},
    {True, True, True, True, True, True, True, True}}

Whehter this should be regarded as a "better" solution of your  
original system or not depends on how you obtained your original  
equations.

Andrzej Kozlowski

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