Iteration of a recurrence solution in R

user_13582844 0
Viewed 166

I'm given a question in R language to find the 30th term of the recurrence relation x(n) = 2*x(n-1) - x(n-2), where x(1) = 0 and x(2) = 1. I know the answer is 29 from mathematical deduction. But as a newbie to R, I'm slightly confused by how to make things work here. The following is my code:

loop <- function(n){
  a <- 0
  b <- 1
  for (i in 1:30){
    a <- b
    b <- 2*b - a
  }
  return(a)
}

loop(30)

I'm returned 1 as a result, which is way off.

In case you're wondering why this looks Python-ish, I've mostly only been exposed to Python programming thus far (I'm new to programming in general). I've tried to check out all the syntax in R, but I suppose my logic is quite fixed by Python. Can someone help me out in this case? In addition, does R have any resources like PythonTutor to help visualise the code execution logic?

Thank you!

recurrence linear-algebra iteration r matrix
4 Answers

I guess what you need might be something like below

loop <- function(n){
  if (n<=2) return(n-1)
  a <- 0
  b <- 1
  for (i in 3:n){
    a_new <- b
    b <- 2*b - a
    a <- a_new
  }
  return(b)
}

then

> loop(30)
[1] 29

If you need a recursion version, below is one realization

loop <- function(n) {
  if (n<=2) return(n-1)
  2*loop(n-1)-loop(n-2)
}

which also gives

> loop(30)
[1] 29
user_121587570

You can solve it another couple of ways.

  1. Solve the linear homogeneous recurrence relation, let

x(n) = r^n

plugging into the recurrence relation, you get the quadratic

r^n-2*r^(n-1)+r^(n-2) = 0

, i.e.,

r^2-2*r+1=0

, i.e.,

r = 1, 1

leading to general solution

x(n) = c1 * 1^n + c2 * n * 1^n = c1 + n * c2

and with x(1) = 0 and x(2) = 1, you get c2 = 1, c1 = -1, s.t.,

x(n) = n - 1

=> x(30) = 29

Hence, R code to compute x(n) as a function of n is trivial, as shown below:

x <- function(n) {
   return (n-1)
}
x(30)
#29

  1. Use matrix powers (first find the following matrix A from the recurrence relation):

    [x(n) x(n-1)]' = [

(The matrix A has algebraic / geometric multiplicity, its corresponding eigenvectors matrix is singular, otherwise you could use spectral decomposition yourself for fast computation of matrix powers, here we shall use the library expm as shown below)

library(expm)

A <- matrix(c(2,1,-1,0), nrow=2)

A %^% 29 %*% c(1,0)  # [x(31) x(30)]T = A^29.[x(2) x(1)]T
#      [,1]
# [1,]   30  # x(31)
# [2,]   29  # x(30)

# compute x(n)
x <- function(n) {
    (A %^% (n-1) %*% c(1,0))[2]
}
x(30)
# 29
user_47061710

You're not using the variable you're iterating on in the loop, so nothing is updating.


loop <- function(n){
  a <- 0
  b <- 1
  for (i in 1:30){
    a <- b
    b <- 2*i - a
  }
  return(a)
}
user_150672970

You could define a recursive function.

f <- function(x, n) {
  n <- 1:n
  r <- function(n) {
    if (length(n) == 2) x[2]
    else r({
      x <<- c(x[2], 2*x[2] - x[1])
      n[-1]
    })
  }
  r(n)
}

x <- c(0, 1)
f(x, 30)
# [1] 29
user_65740380
Related