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