# Função a ser maximizada
f<-function(s){ 
x=s[1];y=s[2]
p=x*exp(-x^2-y^2);p
}

#  Para fazer grafico a função não pode ser definida como vetorial
x <- seq(-2, 2, length = 40)  
y <- seq(-2,2, length = 40)  
g <- function(x,y) {x*exp(-x^2-y^2)} # adptação de f para gerar gráfico
z <- outer(x, y, g) 
persp(x, y, z,theta = 30, phi = 40, zlab = "f(x,y)")->res 

# Derivada da função ou vetor gradiente
gradf<-function(s){ 
x=s[1];y=s[2]
p=exp(-x^2-y^2)*c(1-2*x^2,-2*x*y)
p
}

# --- Método de Euler para resolver a equação u'(t)= grad f(u(t))

b=100
n=10000
h=b/n

u=matrix(0,2,n+1)
u[,1]=c(2,-2)

for (i in 1:n){
u[,i+1]=u[,i]+h*gradf(u[,i])
}

print("Pondo de maximo aproximado"); u[,n+1]

print("Valor maximo aproximado"); f(u[,n+1])

print("Gradiente aproximado no ponto de maximo"); gradf(u[,n+1])

# --- Método de Euler melhorado para resolver a equação u'(t)= -grad f(u(t))

b=100
n=10000
h=b/n

w=matrix(0,2,n+1)
w[,1]=c(-2,2)

for (i in 1:n){
w[,i+1]=w[,i]-h*gradf(w[,i])
}

print("Pondo de minimo aproximado"); w[,n+1]

print("Valor minimo aproximado"); f(w[,n+1])

print("Gradiente aproximado no ponto de minimo"); gradf(w[,n+1])

require(grDevices) 
lines(trans3d(u[1,],u[2,],g(u[1,],u[2,]), res), col = "blue", lwd = 2) # Inclui curva no gráfico 3d.
persp(x, y, z,theta = -60, phi = 60, zlab = "f(x,y)")->res2 
lines(trans3d(w[1,],w[2,],g(w[1,],w[2,]), res2), col = "red", lwd = 2) # Inclui curva no gráfico 3d.
plot(u[1,],u[2,],xlab="u_1(t)",ylab="u_2(t)",col="blue",'l')
plot(w[1,],w[2,],xlab="w_1(t)",ylab="w_2(t)",col="red",'l')



t=seq(0,b,by=h)
plot(t,g(u[1,],u[2,]),ylab="f(u(t))",xlab="t",col="blue",'l')

plot(t,g(w[1,],w[2,]),ylab="f(w(t))",xlab="t",col="red",'l')