Повышение мощности на матрицах

Question

Я пытаюсь создать функцию, которая даст мне значение матрицы, когда она будет поднята до мощности. Это то, что я сделал до сих пор:

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

power <- function(A,n){ 
+ if(n == 0){ 
+ return(diag(4)) 
+ }else{ 
+ return(A%*%A^(n-1)) 
+ } 
+ } 

РЕЗУЛЬТАТ:

> power(A,4) 
[,1] [,2] [,3] [,4] 
[1,] 1 0 0 0 
[2,] 10 1 0 0 
[3,] 46 10 1 0 
[4,] 146 46 10 1 

Это дает другое значение от того, что получает мой калькулятор, и я пытаюсь понять, что я делаю неправильно , Любая помощь приветствуется!

Answer 1

Я предполагаю, что вы хотите сделать muliplication в matrix.You должны сначала сделать умножить матрицу, а затем попытаться умножить их, используя те же полномочия, как вы хотите, так что вы можете сделать две вещи

1. Введите код, чтобы умножить матрицу.
2. Зациклируйте код для умножения.

Answer 2

У вас возникла проблема с тем, как вы вычисляете свой матричный продукт. Вместо этого я использую цикл while внутри вашей функции power(). Он просто умножает матрицу ввода на себя n раз, а затем возвращает результат. Вот базовое R-решение, которое является продолжением направления, в котором вы уже шли.

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

power <- function(A,n){ 
    B <- diag(nrow(A)) 
    if (n == 0) { 
    return(diag(nrow(A))) 
    } else { 
     while (n > 0) { 
      B <- A%*%B 
      n <- n - 1 
     } 
     return(B) 
    } 
} 

> power(A, 4) 
    [,1] [,2] [,3] [,4] 
[1,] 1 0 0 0 
[2,] 8 1 0 0 
[3,] 36 8 1 0 
[4,] 120 36 8 1 

Answer 3

Мы могли использовать %^% из library(expm)

library(expm) 
A%*%(A%^%3) 

Используя это в функции

power <- function(A,n){ 
    if(n == 0){ 
    return(diag(4)) 
    }else{ 
return(A%*%(A%^%(n-1))) 
} 
} 

power(A,4) 
#  [,1] [,2] [,3] [,4] 
#[1,] 1 0 0 0 
#[2,] 8 1 0 0 
#[3,] 36 8 1 0 
#[4,] 120 36 8 1 

Согласно описанию в ?matpow

Вычислить k-ю степень матрицы. В то время как «x^k» вычисляет Элемент wise полномочий, 'x% ^% k' соответствует матрице k - 1 умножения, 'x% % x%% ...% *% x'.

---

Или base R вариант Reduce с %*% (но это будет медленным по сравнению с %^%.

Reduce(`%*%`,replicate(4, A, simplify=FALSE)) 

В функции

power1 <- function(A,n){ 
    if(n == 0){ 
    return(diag(4)) 
    }else{ 
    Reduce(`%*%`,replicate(n, A, simplify=FALSE)) 
    } 
} 

power1(A,4) 
#  [,1] [,2] [,3] [,4] 
#[1,] 1 0 0 0 
#[2,] 8 1 0 0 
#[3,] 36 8 1 0 
#[4,] 120 36 8 1