Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

answered 4 years ago

FrédéricC gravatar image

Use exact fields and polynomial rings, avoid the symbolic ring and numerical reals::

sage: K.<a> = QuadraticField(2)                                                   
sage: x = polygen(K, 'x')                                                          
sage: def yx(x):  
....:     cst = a - 1 
....:     if x == 1: 
....:         return -Infinity 
....:     return (cst * x - 1) / (x - 1) 
....: tmp = x  
....: for t in range(1,9):  
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 ((a - 1)*x - 1)/(x - 1)
2 (a*x - 1)/x
3 (-x - 1)/(x - a - 1)
4 (1/2*a*x - 1)/(x - 1/2*a)
5 ((a + 1)*x - 1)/(x + 1)
6 (-1)/(x - a)
7 (x - 1)/(x - a + 1)
8 x
click to hide/show revision 2
No.2 Revision

Use exact fields and polynomial rings, avoid the symbolic ring and numerical reals::reals:

sage: K.<a> = QuadraticField(2)                                                   
sage: x = polygen(K, 'x')                                                          
sage: def yx(x):  
....:     cst = a - 1 
....:     if x == 1: 
....:         return -Infinity 
....:     return (cst * x - 1) / (x - 1) 
....: tmp = x  
....: for t in range(1,9):  
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 ((a - 1)*x - 1)/(x - 1)
2 (a*x - 1)/x
3 (-x - 1)/(x - a - 1)
4 (1/2*a*x - 1)/(x - 1/2*a)
5 ((a + 1)*x - 1)/(x + 1)
6 (-1)/(x - a)
7 (x - 1)/(x - a + 1)
8 x
click to hide/show revision 3
No.3 Revision

Use exact fields and polynomial rings, avoid the symbolic ring and numerical reals:

sage: K.<a> = QuadraticField(2)                                                   
sage: x = polygen(K, 'x')                                                          
sage: def yx(x):  
....:     cst = a - 1 
....:     if x == 1: 
....:         return -Infinity 
....:     return (cst * x - 1) / (x - 1) 
....: tmp = x  
....: for t in range(1,9):  
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 ((a - 1)*x - 1)/(x - 1)
2 (a*x - 1)/x
3 (-x - 1)/(x - a - 1)
4 (1/2*a*x - 1)/(x - 1/2*a)
5 ((a + 1)*x - 1)/(x + 1)
6 (-1)/(x - a)
7 (x - 1)/(x - a + 1)
8 x

EDIT: and when using this with numerical input, the answers are more reasonable

sage: tmp = 1.5                                                                 
sage: for t in range(1,9): 
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 -0.757359312880715
2 0.747546895706428
3 2.73459080339014
4 0.0765048437046765
5 1.04852813742386
6 -11.6568542494924
7 0.460495713220364
8 1.50000000000000

The correct answer differs at the 6th value from the wrong computation. Not clear how and why.

click to hide/show revision 4
No.4 Revision

Use exact fields and polynomial rings, avoid the symbolic ring and numerical reals:

sage: K.<a> = QuadraticField(2)                                                   
sage: x = polygen(K, 'x')                                                          
sage: def yx(x):  
....:     cst = a - 1 
....:     if x == 1: 
....:         return -Infinity 
....:     return (cst * x - 1) / (x - 1) 
....: tmp = x  
....: for t in range(1,9):  
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 ((a - 1)*x - 1)/(x - 1)
2 (a*x - 1)/x
3 (-x - 1)/(x - a - 1)
4 (1/2*a*x - 1)/(x - 1/2*a)
5 ((a + 1)*x - 1)/(x + 1)
6 (-1)/(x - a)
7 (x - 1)/(x - a + 1)
8 x

EDIT: and when using this with numerical input, the answers are more reasonable

sage: tmp = 1.5                                                                 
sage: for t in range(1,9): 
....:     tmp = yx(tmp) 
....:     print(t, tmp) 
....:                                                                           
1 -0.757359312880715
2 0.747546895706428
3 2.73459080339014
4 0.0765048437046765
5 1.04852813742386
6 -11.6568542494924
7 0.460495713220364
8 1.50000000000000

The correct answer differs at the 6th value from the wrong computation. Not clear how and why.

EDIT:

Comparing the correct answer and the procedure using the symbolic sqrt(2), one can see precision loss happening, with large exponents for the floating point coefficients of sqrt(2)

(-0.757359312880715, 3.0*sqrt(2) - 5.0)
(0.747546895706428, (-8.0*sqrt(2) + 10.0)/(3.0*sqrt(2) - 6.0))
(2.73459080339014, (-150.0*sqrt(2) + 210.0)/(114.0*sqrt(2) - 162.0))
(0.0765048437046765, (-79524.0*sqrt(2) + 112464.0)/(85176.0*sqrt(2) - 120456.0))
(1.04852813742386, (-25732492128.0*sqrt(2) + 36391239360.0)/(39678297120.0*sqrt(2) - 56113585920.0))
(-11.6568542494924, (-(2.51898761979489e+21)*sqrt(2) + 3.56238645536385e+21)/((7.34086788505619e+21)*sqrt(2) - 1.03815549226356e+22))
(0.460495713220364, ((2.61510083243401e+43)*sqrt(2) - 3.69831106420135e+43)/((2.04721262905725e+44)*sqrt(2) - 2.89519586507424e+44))
(1.50000000000000, ((1.55098758806785e+89)*sqrt(2) - 2.19342768211789e+89)/((1.03399172537856e+89)*sqrt(2) - 1.46228512141192e+89))

So this problem is essentially caused by using a mixture of symbolic sqrt(2) and floating points.