Good luck my friend. You are in the undiscovered territory at least as far as I know. If I ever can understand why that is working I will pass it on.
]]>It will be interesting to see if this works for other Taylor expansions (e.g., sin(x), exp(x), etc.).
]]>Have you tried plugging in 1 / 8 into what you are getting and into the Taylor series of ln(1 + x )?
]]>a[0] = (1/8) * ( 1 - (1/8) * ( (1/2) - (1/8) * ( (1/3) - a[3] ))))))...)))))
The coefficient of the "1" is 1/8.
The coefficient of the "1/2" is -(1/8)^2
The coefficient of the "1/3" is (1/8)^3
The Taylor expansion of ln(1+x) = x - (1/2)*x^2 + (1/3)*x^3...
So when x=1/8, for a[0] these 2 series are the same.
]]>In essence, the iterations generate the Taylor series expansion for ln(1+(1/8)). Is this correct?
That I cannot say. As for the 1 / n term, when I first saw this problem I had no inkling it could be part of a larger problem. Recently I found the excellent paper "An example of Round Off Error," by Bruno Guerrieri. Like him I first saw the problem in Cleve Moler's book. Only Bruno went much further with it. Apparently it comes from the attempt to integrate:
The recurrence arises naturally from the attempts to get an accurate answer! That is where the 1 / n comes from. I had to spend a lot of time in remembering and finding that particular recurrence. You see modern compilers and math packages have mechanisms to prevent errors due to smearing. To illustrate the point I had to find one that I was sure would defeat Borland's excellent compiler!
I assume you are using the Taylor expansion of ln(1 + x)?
]]>In essence, the iterations generate the Taylor series expansion for ln(1+(1/8)). Is this correct? If so, this explains why you have the 1/n term. For exp(x), that term would have to be 1/n!, and you shouldn't need the minus sign....?
]]>Remember you start off like this. In this case e is large.
n + e
e/8, e/64, e/512, e/4096...
Pretty obvious that the error is chopped up. Small error = accurate answer!
Remember in the forward direction you needed 32 digits for ln(9/8) and for each iteration, to get 5 or 6 that are accurate. You lost 25 digits. Only natural to say backwards you pick up 25 digits of accuracy. You even can start with a negative number! As long as you iterate enough, any amount of precision is possible in this direction for this problem.
I think this idea is first attributed to Wilkinson, a really brilliant guy. It does seem like magic the first time you see it! It is part of the beauty of numerical analysis!
]]>Although I can believe that errors will grow smaller and smaller when the coefficient is less than 1, it isn't obvious to me why the sequence values are asserting themselves so strongly. But clearly they do...
a[ 51] = 1.00000e+03 (1000)
a[50] = -1.24998e+02 (-124.998)
a[49] = 1.56272e+01 (15.6272)
a[48] = -1.95085e+00 (-1.95085)
a[47] = 2.46460e-01 (0.24646)
a[46] = -2.81479e-02 (-0.0281479)
a[45] = 6.23588e-03 (0.00623588)
a[44] = 1.99829e-03 (0.00199829)
a[43] = 2.59112e-03 (0.00259112)
a[42] = 2.58309e-03 (0.00258309)
a[41] = 2.65330e-03 (0.0026533)
a[40] = 2.71712e-03 (0.00271712)
a[39] = 2.78536e-03 (0.00278536)
a[38] = 2.85696e-03 (0.00285696)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)
a[ 51] = 5.55500e+03 (5555)
a[50] = -6.94373e+02 (-694.373)
a[49] = 8.67991e+01 (86.7991)
a[48] = -1.08473e+01 (-10.8473)
a[47] = 1.35852e+00 (1.35852)
a[46] = -1.67156e-01 (-0.167156)
a[45] = 2.36118e-02 (0.0236118)
a[44] = -1.73701e-04 (-0.000173701)
a[43] = 2.86262e-03 (0.00286262)
a[42] = 2.54915e-03 (0.00254915)
a[41] = 2.65755e-03 (0.00265755)
a[40] = 2.71659e-03 (0.00271659)
a[39] = 2.78543e-03 (0.00278543)
a[38] = 2.85695e-03 (0.00285695)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)
a[ 71] = 5.55500e+03 (5555)
a[70] = -6.94373e+02 (-694.373)
a[69] = 8.67984e+01 (86.7984)
a[68] = -1.08480e+01 (-10.848)
a[67] = 1.35784e+00 (1.35784)
a[66] = -1.67864e-01 (-0.167864)
a[65] = 2.28769e-02 (0.0228769)
a[64] = -9.36541e-04 (-0.000936541)
a[63] = 2.07019e-03 (0.00207019)
a[62] = 1.72535e-03 (0.00172535)
a[61] = 1.80046e-03 (0.00180046)
a[60] = 1.82412e-03 (0.00182412)
a[59] = 1.85532e-03 (0.00185532)
a[58] = 1.88673e-03 (0.00188673)
a[57] = 1.91933e-03 (0.00191933)
a[56] = 1.95307e-03 (0.00195307)
a[55] = 1.98801e-03 (0.00198801)
a[54] = 2.02423e-03 (0.00202423)
a[53] = 2.06179e-03 (0.00206179)
a[52] = 2.10077e-03 (0.00210077)
a[51] = 2.14125e-03 (0.00214125)
a[50] = 2.18332e-03 (0.00218332)
a[49] = 2.22708e-03 (0.00222708)
a[48] = 2.27263e-03 (0.00227263)
a[47] = 2.32009e-03 (0.00232009)
a[46] = 2.36956e-03 (0.00236956)
a[45] = 2.42120e-03 (0.0024212)
a[44] = 2.47513e-03 (0.00247513)
a[43] = 2.53152e-03 (0.00253152)
a[42] = 2.59054e-03 (0.00259054)
a[41] = 2.65237e-03 (0.00265237)
a[40] = 2.71723e-03 (0.00271723)
a[39] = 2.78535e-03 (0.00278535)
a[38] = 2.85696e-03 (0.00285696)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)
a[ 71] = 1.00000e+00 (1)
a[70] = -1.23239e-01 (-0.123239)
a[69] = 1.71906e-02 (0.0171906)
a[68] = -3.37236e-04 (-0.000337236)
a[67] = 1.88039e-03 (0.00188039)
a[66] = 1.63062e-03 (0.00163062)
a[65] = 1.69011e-03 (0.00169011)
a[64] = 1.71181e-03 (0.00171181)
a[63] = 1.73915e-03 (0.00173915)
a[62] = 1.76673e-03 (0.00176673)
a[61] = 1.79529e-03 (0.00179529)
a[60] = 1.82477e-03 (0.00182477)
a[59] = 1.85524e-03 (0.00185524)
a[58] = 1.88674e-03 (0.00188674)
a[57] = 1.91933e-03 (0.00191933)
a[56] = 1.95307e-03 (0.00195307)
a[55] = 1.98801e-03 (0.00198801)
a[54] = 2.02423e-03 (0.00202423)
a[53] = 2.06179e-03 (0.00206179)
a[52] = 2.10077e-03 (0.00210077)
a[51] = 2.14125e-03 (0.00214125)
a[50] = 2.18332e-03 (0.00218332)
a[49] = 2.22708e-03 (0.00222708)
a[48] = 2.27263e-03 (0.00227263)
a[47] = 2.32009e-03 (0.00232009)
a[46] = 2.36956e-03 (0.00236956)
a[45] = 2.42120e-03 (0.0024212)
a[44] = 2.47513e-03 (0.00247513)
a[43] = 2.53152e-03 (0.00253152)
a[42] = 2.59054e-03 (0.00259054)
a[41] = 2.65237e-03 (0.00265237)
a[40] = 2.71723e-03 (0.00271723)
a[39] = 2.78535e-03 (0.00278535)
a[38] = 2.85696e-03 (0.00285696)
a[37] = 2.93235e-03 (0.00293235)
a[36] = 3.01183e-03 (0.00301183)
a[35] = 3.09574e-03 (0.00309574)
a[34] = 3.18446e-03 (0.00318446)
a[33] = 3.27841e-03 (0.00327841)
a[32] = 3.37808e-03 (0.00337808)
a[31] = 3.48399e-03 (0.00348399)
a[30] = 3.59676e-03 (0.00359676)
a[29] = 3.71707e-03 (0.00371707)
a[28] = 3.84571e-03 (0.00384571)
a[27] = 3.98357e-03 (0.00398357)
a[26] = 4.13168e-03 (0.00413168)
a[25] = 4.29123e-03 (0.00429123)
a[24] = 4.46360e-03 (0.0044636)
a[23] = 4.65038e-03 (0.00465038)
a[22] = 4.85348e-03 (0.00485348)
a[21] = 5.07513e-03 (0.00507513)
a[20] = 5.31799e-03 (0.00531799)
a[19] = 5.58525e-03 (0.00558525)
a[18] = 5.88079e-03 (0.00588079)
a[17] = 6.20935e-03 (0.00620935)
a[16] = 6.57677e-03 (0.00657677)
a[15] = 6.99040e-03 (0.0069904)
a[14] = 7.45953e-03 (0.00745953)
a[13] = 7.99613e-03 (0.00799613)
a[12] = 8.61587e-03 (0.00861587)
a[11] = 9.33968e-03 (0.00933968)
a[10] = 1.01962e-02 (0.0101962)
a[ 9] = 1.12255e-02 (0.0112255)
a[ 8] = 1.24857e-02 (0.0124857)
a[ 7] = 1.40643e-02 (0.0140643)
a[ 6] = 1.60991e-02 (0.0160991)
a[ 5] = 1.88209e-02 (0.0188209)
a[ 4] = 2.26474e-02 (0.0226474)
a[ 3] = 2.84191e-02 (0.0284191)
a[ 2] = 3.81143e-02 (0.0381143)
a[ 1] = 5.77357e-02 (0.0577357)
a[ 0] = 1.17783e-01 (0.117783)
Run it in double precision and you will get all the digits!
a[25]=.004 291 231 914 255 080 933 09 7532 748...
]]>However we've lost a[0] - the starting point. It's not obvious to me where you start. If you assume the existence of a fixed point (for n going to infinity), you can solve for the fixed point, but I don't see how to get back down to finite n.
So what do you do?
]]>If multiplying by some constant is unstable when that constant is greater than 1 than dividing by it must be stable.
We can run the recurrence backwards, first you solve for a[n-1].
Now we have the 8 on the bottom as a divisor. The error is going to work for us! What does this difference equation say. Well it runs from the higher numbers downwards ( backwards direction ). We need to start from say a[50] then we would get a[49], a[48], a[47]... all the way down to a[25].
Only one problem, how do we know what a[50] is? We need a starting value as the initial condition. How do we get that when we could not even get a[25]?
What would you do?
Yes , that is almost right. Not only transcendentals and irrational numbers but remember post #28? Good old fractions like 1 / 3 and 1 / 10 ... are not representable in a finite number of bits.
Anytime a multiplication occurs 101 * .123456 = 12.4690, actual answer 12.469006. We lost 2 digits, they were smeared out.
Anytime you have a number n, unless you are dealing with integers you have a small error called little e. You do not really have n but instead n + e. Everytime you do a calculation the little e is either no trouble because it remains small or...
Here is what happened.
1st iteration) 8 (n + e) = 8n + 8e. For the sake of brevity I will just chart the e.
2nd iteration) 8 (n + 8 e) = 64e
3rd iteration) 8 (n + 64 e) = 512e
3th iteration) 8 (n + 512 e) = 4096e etc.
The small e is now 4000 times larger than it was. It started at e = .0000005 and now it is approximately .002. In other words it is affecting the second and third digits after the decimal point. You already have lost 4 - 5 significant digits!
That is the reason the forward direction is impossible. Since you lost 1 digit for every iteration you would have had to start with ln(9/8) to 31-32 digits and all subsequent calculations to 31-32 digits to have 5 or 6 correct at the end.
Do you follow me so far?
]]>Okay, notice that the numbers are oscillating plus and minus. The actual answer is:
a[25]=.004 291 231 914 255 080 933 09 7532 748... to 30 digits.
The double precision answer is - 45312.4, 10 million times too large and having the wrong sign. What went wrong. Nothing, no matter what the precision the round off error will eventually catch up to you. The problem is that the difference equation, is ill conditioned or pathogenic.
This can be added to the above 6 in post #28. 7) Computers cannot even input numbers correctly.
Remember a[0] = ln(9 / 8) you put it in, in double precision ( around 17 digits ). But that was not really ln(9/8), it had a small error. Each time time you iterated throught the recurrence you magnified that error by 8 until the error swamped out the correct digits, that is what smearing is.
That is why there are no correct digits in your answer. In order to get 1 correct digit you would have needed to work with 25 digits of precision.
Do you understand why the 8 as a multiplier is causing the problem. If you do not understand tell me now or otherwise the fix will have no meaning to you.
]]>