Recurrence Relation for a Decreasing Exponential Has Increasing Values

Question

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.

Answer 1

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).