Power Method Animation

Code for drawing

power <- function(r, x1 = c(1, 0), x2 = c(1, 0), eps = 0.001, itmax = 100) {
    m <- matrix(c(1, r, r, 1), 2, 2)
    n <- solve(m)
    itel <- 1
    repeat {
        y <- colSums(x1 * m)
        z <- colSums(x2 * n)
        plot(ellipse(matrix(c(1, 0, 0, 1), 2, 2), t = 1), xlim = c(-2, 2), ylim = c(-2, 
            2), type = "l", asp = 1, lwd = 3, col = "BLUE")
        lines(matrix(c(0, x1[1], 0, x1[2]), 2, 2), lwd = 2)
        lines(matrix(c(0, x2[1], 0, x2[2]), 2, 2), lwd = 2)
        arrows(0, 0, y[1], y[2], col = "RED", lwd = 2)
        arrows(0, 0, z[1], z[2], col = "GREEN", lwd = 2)
        s <- sqrt(sum(y^2))
        u <- sqrt(sum(z^2))
        text(y[1] + 0.1, y[2] + 0.1, format(s, digits = 2))
        text(z[1] - 0.1, z[2] - 0.1, format(1/u, digits = 2))
        y <- y/s
        z <- z/u
        ops <- max(max(abs(x1 - y)), max(abs(x2 - z)))
        if ((itel == itmax) || (ops < eps)) {
            (break)()
        }
        x1 <- y
        x2 <- z
        itel <- itel + 1
    }
}

Make the animation

library(ellipse)
library(animation)
r <- 0.2
ani.options(loop = FALSE)
saveHTML(power(r), interval = 0.8)

## HTML file created at: index.html

Show the animation

system("open index.html")