Re: A problem about numerical precision

*To*: mathgroup at smc.vnet.net*Subject*: [mg52646] Re: [mg52609] A problem about numerical precision*From*: DrBob <drbob at bigfoot.com>*Date*: Sun, 5 Dec 2004 02:08:53 -0500 (EST)*References*: <200412040907.EAA13394@smc.vnet.net>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

This is a situation where Mathematica's precision arithmetic seems to get us into trouble. Preliminaries: short[n_:7][expr_] := Short[expr, n] delta = 1/100; a = 9/10; exactB = -3*Sqrt[2]/10; b = exactB; A = {{1, 0}, {b, a}}; B = {{-1, 1}, {-b, b}}; f[x_, y_] := If[Abs[y - x] > delta, If[ones == 0, Print@{zeroes, Abs[x - y]}]; ones++; 1, zeroes++; 0]; Now consider the eigenvalues of A and B: Eigenvalues /@ {A, B} {{1, 9/10}, {(1/10)*(-10 - 3*Sqrt[2]), 0}} Depending on the starting point, the series can converge on {0,0} or any multiple of A's first eigenvector, so long as the coordinates of the vector are no more than delta apart. Otherwise B.{x,y} is added, which causes the series to keep moving around rather than settling on an eigenvector of A. Eigenvectors[A] {{-(1/(3*Sqrt[2])), 1}, {0, 1}} Notice that if x starts out at infinite precision (as it does), it will remain at infinite precision until B.{x,y} is added, at which point we have a situation where the new x is the old y. A . {x, y} B . {x, y} % + %% {x, -((3*x)/(5*Sqrt[2])) + (9*y)/10} {-x + y, (3*x)/(5*Sqrt[2]) - (3*y)/(5*Sqrt[2])} {y, (9*y)/10 - (3*y)/(5*Sqrt[2])} But that happened by computing newY = (old x) + (oldY - oldX) = oldY, which may be a numerically unstable computation. The test in f is inherently vulnerable to numerical instability as well, and B has an eigenvalue less than -1, so much depends on how often the test evaluates to True. Abs[x-y] can be very close to delta on either side of it, so this is sensitive to starting points and the precision of the inputs. We can't evaluate 3000 terms with b in its exact state, so we approximate it. For instance, b = N@exactB; ones = count = 0; {it, {zeroMachine, oneMachine}} = Reap at NestList[A.# + f @@ #B.# &, -delta/2{1, a + b}, 3000 - 1]; {count, ones} (temp = Take[it, -1000]) // short[] Map[Precision, it, {2}] // short[] ListPlot[temp, PlotRange -> All]; {2999, 644} {{0.00660334,-0.00277839},{0.00660334,-0.00530212},{-0.00530212,-0.00252241},\ \[LeftSkeleton]995\[RightSkeleton],{-0.00535461,-0.00254738},{-0.00535461,-0.\ 0000208737}} {{â??,MachinePrecision},{MachinePrecision,MachinePrecision},{MachinePrecision,\ MachinePrecision},{MachinePrecision,MachinePrecision},\[LeftSkeleton]2993\ \[RightSkeleton],{MachinePrecision,MachinePrecision},{MachinePrecision,\ MachinePrecision},{MachinePrecision,MachinePrecision}} B.{x,y} was added 644 times, and we didn't get convergence in 3000 iterations. If we increase the precision of b, we get very different results: b = N[exactB, 200]; ones = count = 0; {it, {zero200, one200}} = Reap@NestList[ A.# + f @@ #B.# &, -delta/2{1, a + b}, 3000 - 1]; {count, ones} (temp = Take[it, -1000]) // short[] Map[Precision, it, {2}] // short[] ListPlot[temp, PlotRange -> All]; {2999,335} {{0.015,-0.065},{0.015,-0.065},\[LeftSkeleton]997\[RightSkeleton],{0.015,-0.\ 065}} {{â??,200.05},{â??,198.009},{â??,200.088},{â??,200.414},\[LeftSkeleton]2992\ \[RightSkeleton],{2.0986,2.0986},{2.0986,2.0986},{2.0986,2.0986},{2.0986,2.\ 0986}} B.{x,y} was added only 335 times, convergence was achieved before the 2000th iteration, and the answers have gradually lost precision. Now, let's check to see when and how the first major difference occurred in the two series: when B.{x,y} was first added in one series but not the other, and vice versa. A little bit of sleuthing gives the answer: first=First@ Complement[ Union[zeroMachine[[All, 1]],zero200[[All,1]]],Intersection[zeroMachine[[All,1]],zero200[[All, 1]]]] Position[zero200,first,Infinity,1] Position[oneMachine,first,Infinity,1] zero200[[1207,2]]//InputForm 1536 {{1207,1}} {{330,1}} 0.0100057828055487347`5.133659435497401 oneMachine[[330,2]]//InputForm 0.010005782805548664 At iteration number 1536, both series have virtually the same Abs[x-y], but in one case it has about 5 digits of precision and in the other, it has $MachinePrecision. That's why B.{x,y} is added in the latter case, but not the former -- because the test comes out False when there isn't enough precision to say that .010000578 > .01. So, arguably, increasing precision caused this computation to go wrong, when lower precision would make it come out right, because increasing the precision of b put us into a different numerical model that let precision atrophy away. Mathematica is needlessly pessimistic about how accurate the result is, at this point. Which plot is more correct, then? Let's set precision at a FIXED level of 200, and see what happens: b = N[exactB, 200]; ones = count = 0; Block[{$MaxPrecision = 200, $MinPrecision = 200}, {it, {zero200, one200}} = \ Reap at NestList[A.# + f @@ #B.# &, -delta/2{1, a + b}, 3000 - 1];] {count, ones} (temp = Take[it, -1000]) // N // short[] Map[Precision, it, {2}] // short[] ListPlot[temp, PlotRange -> All]; {2999,644} {{0.00660334369961137174735312435438344177537767846967480208854829063743531877\ 669654426058000132359435471231012469047871989381708367247411510957737239871571\ 86279055630437225780845060015289235201211328100414,-\[LeftSkeleton]227\ \[RightSkeleton]},\[LeftSkeleton]998\[RightSkeleton],{-\[LeftSkeleton]218\ \[RightSkeleton],-\[LeftSkeleton]228\[RightSkeleton]}} {{â??,200.},{200.,200.},{200.,200.},\[LeftSkeleton]2994\[RightSkeleton],{200.,\ 200.},{200.,200.},{200.,200.}} We're back to the original plot, almost exactly! Bobby On Sat, 4 Dec 2004 04:07:48 -0500 (EST), Guofeng Zhang <guofengzhang at gmail.com> 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 program. 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 > > > > -- DrBob at bigfoot.com www.eclecticdreams.net

**References**:**A problem about numerical precision***From:*Guofeng Zhang <guofengzhang@gmail.com>