Re: Re: A problem of numerical precision
- To: mathgroup at smc.vnet.net
- Subject: [mg52651] Re: [mg52642] Re: [mg52610] A problem of numerical precision
- From: DrBob <drbob at bigfoot.com>
- Date: Tue, 7 Dec 2004 04:09:36 -0500 (EST)
- References: <200412030815.DAA25377@smc.vnet.net> <200412040907.EAA13400@smc.vnet.net> <200412050708.CAA26788@smc.vnet.net> <opsijs5ropiz9bcq@monster.cox-internet.com>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
Oops, I was using an exact A and B also. Correcting that changes details, but the 200-digit fixed-precision result still matches the machine-precision result -- in both cases the test in f oscillates, and the result of f is precisely the same for any iteration.
But f gives different results when we start at 200 digits and let precision atrophy.
Bobby
On Sun, 05 Dec 2004 12:47:41 -0600, DrBob <drbob at bigfoot.com> wrote:
> You used exact versions of A and B in all your trials, where the OP (I think) and I (definitely) defined A and B based on the approximation to b.
>
> I believe Peter Pein did the opposite thing -- he used the lower-precision version of A and B when using the higher-precision b to define it[[1]].
>
> So we all got different results, just as we could expect to do.
>
> Short of using exact arithmetic throughout, I think the fixed-precision method is authoritative. The result was indistinguishable from the machine-precision values, for a very good reason -- errors don't actually accumulate in this calculation, as Mathematica's precision arithmetic assumes they do. Using fixed-precision OR machine-precision sidesteps that arithmetic and prevents the test in f from going wrong.
>
> If we reach an eigenvector corresponding to A's first eigenvalue, we should STAY there (unless the test in f comes True). If we're anywhere else, A.{x,y} tends to diminish error, not increase it. Even B.{x,y} decreases error for many inputs, and the term isn't usually added anyway. If it isn't added when it shouldn't be, no problem!
>
> Hence numerical instability comes almost exclusively from the test in f, and that's affected far less by accumulated error than by artificially diminished precision.
>
> Bobby
>
> On Sun, 5 Dec 2004 02:08:38 -0500 (EST), Daniel Lichtblau <danl at wolfram.com> wrote:
>
>> Guofeng Zhang wrote:
>>>
>>> Hi,
>>>
>>> I met one problem when I do some iteration: By increasing numerical
>>> precision, the results are so different from the original one! I
>>> don't know what went wrong, and hope to get some answers.
>>>
>>> The code is
>>>
>>> delta = 1/100;
>>> a = 9/10;
>>> b = -3*1.4142135623730950/10;
>>> A = { {1,0}, {b,a} };
>>> B= { {-1,1}, {-b,b} };
>>>
>>> f[v_,x_] := If[ Abs[ (10^20)*v-(10^20)*x]>(10^20)*delta, 1, 0 ];
>>>
>>> M = 3000;
>>> it = Table[0, {i,M}, {j,2} ];
>>> it[ [1] ] = { -delta/2, (a+b)*(-delta/2) };
>>>
>>> For[ i=1, i<M, it[ [i+1]] = A.it[[i]]+f[ it[ [i,1] ], it[ [i,2] ]
>>> ]*B.it[ [i] ]; i++ ];
>>> temp = Take[it, -1000];
>>>
>>> ListPlot[ temp ];
>>>
>>> By evaluating this, I got an oscillating orbit. I got the same using
>>> another system. However, if I increase the numerical precision by using
>>> b = -3*1.4142135623730950``200/10;
>>> to substitute the original b, the trajectory would converge to a fixed
>>> point instead of wandering around!
>>>
>>> I don't know why. I am hoping for suggestions. Thanks a lot.
>>>
>>> Guofeng
>>>
>>
>>
>> You are seeing the effect of plotting "zeros" of negative accuracy.
>>
>> Try it as below.
>>
>> delta = 1/100;
>> a = 9/10;
>> b = -3/10*Sqrt[2];
>> A = {{1,0}, {b,a}};
>> B = {{-1,1}, {-b,b}};
>> f2[v_,x_] := UnitStep[Abs[v-x]-delta]
>> M = 3000;
>>
>> itm = Table[0, {i,M}, {j,2}];
>> itm[[1]] = N[{-delta/2, (a+b)*(-delta/2)}, MachinePrecision];
>> Do[itm[[i+1]] = A.itm[[i]] + f2[itm[[i,1]],itm[[i,2]]]*B.itm[[i]],
>> {i,1,M-1}]
>>
>> it200 = Table[0, {i,M}, {j,2}];
>> it200[[1]] = N[{-delta/2, (a+b)*(-delta/2)}, 200];
>> Do[it200[[i+1]] = A.it200[[i]] + f2[it200[[i,1]],it200[[i,2]]]*B.it200[[i]],
>> {i,1,M-1}]
>>
>> it500 = Table[0, {i,M}, {j,2}];
>> it500[[1]] = N[{-delta/2, (a+b)*(-delta/2)}, 500];
>> Do[it500[[i+1]] = A.it500[[i]] + f2[it500[[i,1]],it500[[i,2]]]*B.it500[[i]],
>> {i,1,M-1}]
>>
>> (These could also be computed with FoldList).
>>
>> Now look at the last elements, or take the Precision, of it500. You will
>> find that fewer than 13 digits are correct. At lower initial precision
>> you get what amount to indeterminate values (except we have some
>> information bounding their magnitudes). Using machine arithmetic there
>> is no way to tell from the result whether it is "correct" or just how
>> far off it might be. It turns out that they are just fine in comparison
>> to it500. So plots of itm and it500 are reliable.
>>
>> The last 1000+ values for it200 are a different matter. They all appear
>> to be zero. For example:
>>
>> In[156]:= it200[[1888]] == {0,0}
>> Out[156]= True
>>
>> A closer inspection reveals that they have negative accuracy. The
>> 1888th, for example, could be anywhere up to around 262 digits in
>> magnitude (and it gets worse as we move further along in it200).
>>
>> In[157]:= InputForm[it200[[1888]]]
>> Out[157]//InputForm= {0``-262.03379085102995, 0``-261.84220176123233}
>>
>> So they plot as zeros, but that's not a very reliable indication of the
>> actual values.
>>
>>
>> Daniel Lichtblau
>> Wolfram Research
>>
>>
>>
>>
>>
>
>
>
--
DrBob at bigfoot.com
www.eclecticdreams.net
- References:
- A problem of numerical precision
- From: Guofeng Zhang <guofengzhang@gmail.com>
- Re: A problem of numerical precision
- From: Daniel Lichtblau <danl@wolfram.com>
- A problem of numerical precision