# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo - Metodo de Newton-Raphson
#==============================================================================
fx <- function(x){
	if(length(x)!=3 || !is.vector(x)){
		print("O vetor deve ter somente 3 elementos!")
		return(NA)
	}	
	y    <- rep(NA,3)
	y[1] <- exp(x[1])-2
	y[2] <- 5*x[3]-4
	y[3] <- 4*x[1]*x[2]-2*x[3]-6	
	return(y)
}

#==============================================================================
# Jacobiano
#==============================================================================
jacobiano.fx <- function(x){
	if(length(x)!=3 || !is.vector(x)){
		print("O vetor deve ter somente 3 elementos!")
		return(NA)
	}	
	J      <- matrix(0,3,3)
	J[1,1] <- exp(x[1])
	J[2,3] <- 5
	J[3,1] <- 4*x[2]	
	J[3,2] <- 4*x[1]	
	J[3,3] <- -2
	return(J)
}

#==============================================================================
# Newton-Raphson
#==============================================================================
newton.raphson.fx <- function(x0=c(1,1,1),nmax.iter=1000,
                                         tolerancia=1.0e-10){
	if(length(x0)!=3 || !is.vector(x0)){
		print("O vetor deve ter somente 3 elementos!")
		return(NA)
	}
	x     <- matrix(NA,(nmax.iter+1),3) 
	x[1,] <- x0
	for(i in 2:(nmax.iter+1)){
		print(paste("Iteracao:",i-1,sep=" "))
		J     <- jacobiano.fx(x[i-1,])
		x[i,] <- x[i-1,] - solve(J)%*%fx(x[i-1,])
		if(sqrt(sum( (x[i,] - x[i-1,])^2 ) ) <tolerancia) 
			return(x[i,])
	}
	print("Numero maximo de iteracoes atingido!")
	return(x[(nmax.iter+1),])
}

#==============================================================================
# Exemplo
#==============================================================================
# Metodo geral com padroes definidos pela funcao 
#==============================================================================
print("===============================")
z1 <- newton.raphson.fx()
print("===============================")
print(z1)
print("===============================")
print(fx(z1))
print("===============================")
print(jacobiano.fx(z1))

[1] "==============================="
[1] "Iteracao: 1"
[1] "Iteracao: 2"
[1] "Iteracao: 3"
[1] "Iteracao: 4"
[1] "Iteracao: 5"
[1] "==============================="
[1] 0,6931472 2,7411206 0,8000000
[1] "==============================="
[1] 0 0 0
[1] "==============================="
         [,1]     [,2] [,3]
[1,]  2,00000 0,000000    0
[2,]  0,00000 0,000000    5
[3,] 10,96448 2,772589   -2

#==============================================================================
# Exemplo
#==============================================================================
# Utilizando diferentes valores iniciais
#==============================================================================
print("===============================")
z2 <- newton.raphson.fx(c(10,10,10))
print("===============================")
z3 <- newton.raphson.fx(c(-1,-1,-1))
print("===============================")
print(z1)
print("===============================")
print(z2)
print("===============================")
print(z3)

[1] "==============================="
[1] "Iteracao: 1"
[1] "Iteracao: 2"
[1] "Iteracao: 3"
[1] "Iteracao: 4"
[1] "Iteracao: 5"
[1] "Iteracao: 6"
[1] "Iteracao: 7"
[1] "Iteracao: 8"
[1] "Iteracao: 9"
[1] "Iteracao: 10"
[1] "Iteracao: 11"
[1] "Iteracao: 12"
[1] "Iteracao: 13"
[1] "Iteracao: 14"
[1] "Iteracao: 15"
[1] "==============================="
[1] "Iteracao: 1"
[1] "Iteracao: 2"
[1] "Iteracao: 3"
[1] "Iteracao: 4"
[1] "Iteracao: 5"
[1] "Iteracao: 6"
[1] "Iteracao: 7"
[1] "Iteracao: 8"
[1] "Iteracao: 9"
[1] "Iteracao: 10"
[1] "==============================="
[1] 0,6931472 2,7411206 0,8000000
[1] "==============================="
[1] 0,6931472 2,7411206 0,8000000
[1] "==============================="
[1] 0,6931472 2,7411206 0,8000000

#==============================================================================
# Nao deixe de verificar as condicoes de convergencia do algoritmo de 
# Newton-Raphson
#==============================================================================

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