# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Funcao para a solucao de sistemas lineares pelo metodo iterativo de Jacobi
#==============================================================================
jacobi <- function(A,b,maxit,x0){
    D         <- as.vector(diag(A))
    LU        <- A    # L + U
    diag(LU)  <- 0    # Zeros na diagonal de L+U
    x         <- x0
    x.new     <- matrix(NA,maxit+1,ncol(A))
    x.new[1,] <- x0
    for(i in 1:maxit){
        x           <- (1/D)*(b-LU%*%x)
        x.new[i+1,] <- x
    }
    return(x.new)
}

#==============================================================================
# Exemplo 
#==============================================================================
L <- matrix(c(1,0,0,2,1,0,2,4,1),3,3,byrow=T)
b <- c(3,4,5)
print("===============================")
print(L)
print("===============================")
print(b)

[1] "==============================="
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    2    1    0
[3,]    2    4    1
[1] "==============================="
[1] 3 4 5

#==============================================================================
# Solucao utilizando a funcao do R
#==============================================================================
result <- solve(L,b)
result  # 3 -2 7

1. 3
2. -2
3. 7

#==============================================================================
# Para diferentes valores iniciais
#==============================================================================
print("===============================")
print(jacobi(L,b,5,c(0,0,0)))
print("===============================")
print(jacobi(L,b,5,c(10,10,10)))

[1] "==============================="
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    3    4    5
[3,]    3   -2  -17
[4,]    3   -2    7
[5,]    3   -2    7
[6,]    3   -2    7
[1] "==============================="
     [,1] [,2] [,3]
[1,]   10   10   10
[2,]    3  -16  -55
[3,]    3   -2   63
[4,]    3   -2    7
[5,]    3   -2    7
[6,]    3   -2    7

# rm(list=ls()) 
#==============================================================================
# Fim
#==============================================================================