n <- 35
F <- rep(0,n)
N <- rep(0,n)
F[1] <- 1
F[2] <- 1/3
for (k in 3:n) F[k] <- (10/3)*F[k-1]- F[k-2]
F
N <- seq(from=1, to=n, by=1)
If you are not familiar with solving linear recurrence equation, it does not matter at all. Anyway, we can get the result of F[n]=3^(1-n) from solving recurrence equation as above, i.e F[n]=(10/3)F[n-1]-F[n-2], f1=1, f2=1/3.
For this reason, by using
plot (N, F,type="l")
we can expect the graph of "3^(1-n)" as known as exponential function.
However, the output is different from the expected. In comparison with the output by
curve(3^(1-x),0,35, add=TRUE, col='blue')
As you know, 3^(1-x) is monotonic decrease function. In spite of expectation, we only get the graph which is increasing in late calculation.
F[18]>F[19]
TRUE
F[19]>F[20]
FALSE
What happened? In common sense, all out put of "F[n]>F[n+1]" should be TRUE.
If I increase the number which is allocated to "n" to 50 from 35,
n <- 50
plot (N, F,type="l")
the shape of graph become totally strange.
I am guessing that the reason is based on the "Double precision binary floating point" (http://en.wikipedia.org/wiki/Double_precision). In my opinion, R allocate the number which is smaller than 0.0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0001 * 2^(-52) (there are 52 zeros) as more large number inversely in recurrence relation.
However, I do not know whether my assumption is true. Even if my assumption is true, why does R allocate very small number as more large number inversely, ONLY in "recurrence relation" not for general function such as 3^(n-1)? Furthermore, in the case of "n=50" why does R change the shape of graph totally?
Could you help me?
Thank you in advance.
This has nothing to do with R, per-se, and everything to do with floating-point values as represented by your computer.
Recurrence relations are like differential equations, and the problem is stated in two parts - the relation and the initial conditions. Change the initial conditions, and you have a different solution.
Note that with initial conditions F[1] <- 1; F[2] <- 3
, the solution is 3^(x-1)
(stated without proof, but it is easy to verify). An increasing exponential function.
Next, note the ratio between elements (it is mildly instructive to also look at intermediate values of H
here):
H <- tail(F, -1) / head(F, -1)
c(head(H, 1), tail(H, 1))
## [1] 0.3333333 3.0000000
You are transitioning between the solution f(x) = 3^(1-x) and f(x) = 3^(x-k) (for some constant k - it is not 1 here, but it is pointless to compute it exactly).
The reason is that when you subtract F[k-2], the arithmetic is not exact, so you do not subtract quite enough at each stage, and it is as if you had a higher initial condition for the exact solution at that stage.
It is valid to give the first N points of F, then use the recurrence relation to solve at that stage. This gives a sequence of functions. And that is what happens when it is computed numerically - it is a different set of initial conditions at each computation.
You're actually computing the solution for f(x) = (10/3)f(x-1) - f(x-2) + e(f(x-2)) where e(x) > 0 for all x (and represents the bits that fall off the end in the subtraction).