Re: Re: Cannot NSolve a system of equations
- To: mathgroup at smc.vnet.net
- Subject: [mg88988] Re: [mg88944] Re: Cannot NSolve a system of equations
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Fri, 23 May 2008 03:06:15 -0400 (EDT)
- References: <g0uahi$4tk$1@smc.vnet.net> <g11rbc$agd$1@smc.vnet.net> <200805220634.CAA22304@smc.vnet.net> <B54B57F9-499C-45E8-AA74-A664AA9EFA83@mimuw.edu.pl>
A correction. The situation looks much brighter. With just two changes in your code: phi = Rationalize[0.74615385, 0]; and zet = 5 taubar/(y1^(1 + phi) + y2^(1 + phi) + y3^(1 + phi) + y4^(1 + phi) + y5^(1 + phi)); I get: sol = N[NSolve[N[eqns1, 500], vars]] {{x -> 0.622656, x1 -> 0.646962, x2 -> 0.649281, x3 -> 0.645755, x4 -> 0.634993, x5 -> 0.536288, abar -> 0.154881, roverw -> 0.212256}, {x -> 1., x1 -> 1., x2 -> 1., x3 -> 1., x4 -> 1., x5 -> 1., abar -> 0.160682, roverw -> 0.}} Both solutions satisfy your requirements but they are different from the one you got. Moreover, if I change phi = Rationalize[0.74615385, 0] to phi = Rationalize[0.7461538, 0]; I get: {{x -> 0.6226558058693201, x1 -> 0.6469616894031311, x2 -> 0.649281289391567, x3 -> 0.6457551146893872, x4 -> 0.6349933630627287, x5 -> 0.5362875727997868, abar -> 0.15488096993390985, roverw -> 0.2122561091985074}, {x -> 1., x1 -> 1., x2 -> 1., x3 -> 1., x4 -> 1., x5 -> 1., abar -> 0.16068192162403339, roverw -> 0.}} {{x -> 0.622656, x1 -> 0.646962, x2 -> 0.649281, x3 -> 0.645755, x4 -> 0.634993, x5 -> 0.536288, abar -> 0.154881, roverw -> 0.212256}, {x -> 1., x1 -> 1., x2 -> 1., x3 -> 1., x4 -> 1., x5 -> 1., abar -> 0.160682, roverw -> 0.}} which is practically identical. So the instability you mentioned has effectively disappeared. Andrzej Kozlowski On 22 May 2008, at 20:56, Andrzej Kozlowski wrote: > 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
- Re: Cannot NSolve a system of equations