Re: Re: Cannot NSolve a system of equations

*To*: mathgroup at smc.vnet.net*Subject*: [mg88995] Re: [mg88944] Re: Cannot NSolve a system of equations*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Fri, 23 May 2008 03:07:32 -0400 (EDT)*References*: <g0uahi$4tk$1@smc.vnet.net> <g11rbc$agd$1@smc.vnet.net> <200805220634.CAA22304@smc.vnet.net>

Note that: phi = Rationalize[0.74615385] so rationalizing phi in this way does nothing at all. Moreover, the fact that your zet is a machine precision number means that the computation is done with machine precision numbers. There is, of course, no unique way to rationalize approximate numbers. When I use Rationalize[ ,0] to force everything in your equations to be rational, I do not get any solutions which satisfy your condition that all the values should lie between 0 and 1. But what is more strange: when simply copy and past your code into a new notebook and evaluate it with fresh kernel I do not get your solution. This may be platform dependent (I am using 6.02 on a MacBook). I think the fact that the computation with machine precision is numerically unstable is not surprising. What is surprising is that the answer that you obtained appears quite accurate - even though it was done essentially with machine precision numbers. That may not be entirely true, since some of the equations involve rational numbers so some of the intermediate computations may well have bee done with high precision, so this could explain why your answer seems quite accurate. Still, the fact that, even running your code without any changes I get a different answer is strange, to say the least. I hope this is verified by other people using Macs. Andrzej Kozlowski PS. When I run your code the only answer that I get that satisfies your conditions that everything lies between 0 and 1 is {x -> 1., x1 -> 1., x2 -> 1., x3 -> 1., x4 -> 1., x5 -> 1., abar -> 0.160682, roverw -> 0.} which is quite far form yours. As I wrote earlier, if the equations are genuinely rationalized, none of the solutions satisfies your conditions. Andrzej Kozlowski On 22 May 2008, at 15:34, murat.koyuncu at gmail.com wrote: > I sincerely thank you all for your help. It cleared up a lot of things > for me. I had tried Rationalize but didn't know about precision > issues. > > Now I have a follow-up question. I still cannot trust my results 100% > because when I slightly change one of my parameters, solution set > changes drastically. My theory says that resulting variables should > between 0 and 1, so I am only interested in that range. And sometimes > I don't even get any such solutions. > > For example, here is my model (very much polished, thanks to your > suggestions): > > Unprotect[In, Out]; Clear[In, Out]; ClearAll["Global`*"]; > phi = Rationalize[0.74615385]; > eta = 7/4; > alpha = 16/25; > taubar = 13/100; > {y1, y2, y3, y4, y5} = > Rationalize[{0.235457064, 0.512465374, 0.781779009, 1.109572176, > 2.360726377}, 0]; > zet = N[5 taubar/(y1^(1 + phi) + y2^(1 + phi) + y3^(1 + phi) + y4^( > 1 + phi) + y5^(1 + phi))]; > tau1 = zet y1^phi; > tau2 = zet y2^phi; > tau3 = zet y3^phi; > tau4 = zet y4^phi; > tau5 = zet y5^phi; > 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)/(5 x), > roverw == (1 - x) (1 - alpha)/alpha}; > vars = {x, x1, x2, x3, x4, x5, abar, roverw}; > sol = NSolve[N[eqns1, 500], vars]; > InputForm[N[sol]] > > One of the solutions I get is between 0 and 1: > {x -> 0.5004885149810059, x1 -> 0.6893393376625371, > x2 -> 0.6337652704581314, x3 -> 0.5731109245369462, > x4 -> 0.4917632705087048, x5 -> 0.1144637717387098, > abar -> 0.1244737777525507, roverw -> 0.2809752103231842} > > So I am happy. But when I drop the last digit of the first parameter, > i.e. phi = Rationalize[0.7461538];, the new set of solutions do not > have any results compatible with my assumptions: > > {{x -> 0.2949617829869077, x1 -> -1.9491786143565832, x2 -> > 10.96762641745164, x3 -> 96.0369893781559, x4 -> -73.18097932930044, > x5 -> -30.399648937016046, abar -> -5.824073507107748, roverw -> > 0.39658399706986447}, > {x -> 2.0689803436369942, x1 -> -8.65276583855146, x2 -> > -18.39904428993495, x3 -> -33.55295774059405, x4 -> > -69.9736576669751, > x5 -> 140.92332725424055, abar -> 2.4458078645586667, roverw -> > -0.6013014432958093}, > {x -> 0.6213605439915423, x1 -> -12.823085779076337, x2 -> > 41.57610980587073, x3 -> -6.337885849417286, x4 -> > -8.665252197121296, > x5 -> -10.64308326029811, abar -> -0.7583206767716181, roverw -> > 0.21298469400475745}, > {x -> 0.6854511826071373 - 0.04362902433398709*I, x1 -> > 0.59964456896311 + 0.08549874839569205*I, > x2 -> 0.6411455792212575 + 0.04079939260558141*I, x3 -> > 0.6760747168284579 - 0.003707834204090894*I, > x4 -> 0.7125762414563436 - 0.05913956075670094*I, x5 -> > 0.7978148065665167 - 0.28159586771041706*I, > abar -> 0.17051392124951886 - 0.016129276403199968*I, roverw -> > 0.17693370978348527 + 0.024541326187867737*I}, > {x -> 0.6854511826071373 + 0.04362902433398709*I, x1 -> > 0.59964456896311 - 0.08549874839569205*I, > x2 -> 0.6411455792212575 - 0.04079939260558141*I, x3 -> > 0.6760747168284579 + 0.003707834204090894*I, > x4 -> 0.7125762414563436 + 0.05913956075670094*I, x5 -> > 0.7978148065665167 + 0.28159586771041706*I, > abar -> 0.17051392124951886 + 0.016129276403199968*I, roverw -> > 0.17693370978348527 - 0.024541326187867737*I}, > {x -> 0.9999999999989753, x1 -> 0.9999999999982316, x2 -> > 0.9999999999985612, x3 -> 0.999999999998809, x4 -> > 0.9999999999991211, > x5 -> 1.0000000000001525, abar -> 0.16068192162412143, roverw -> > 0.}} > > And this is just one of the instances that this happens. I solve the > same set of equations for many parameter sets. Sometimes I get a > solution in the [0,1] range, but it is so different than the previous > one, it doesn't make sense. Or sometimes, I have to tweak parameter > values (drop a digit here, add another one there) to get a meaningful > solution. I need a set of results over different parameter values that > I can compare and contrast, but it is hard to do that when my results > suddenly disappear or jump to an improbable value. > > Now I know that this is probably the most Mathematica can do for me, > but is there a way to make my system more "stable"? Maybe play with my > equations a little bit? Or put another way, what makes my system so > unstable? > > Again, any suggestions will be very much appreciated. Thanks! >

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