# Minimizando funções 


mf<-function(s,r){z=sqrt(s*r)*exp(-sqrt(s^3)/2-r^2)# f(x,y) para ser maximizada.
z
}

require(grDevices) # Gráfico de f(x,y)
x <- seq(0,2.5, by=0.05)
y <- seq(0,1,by=0.05)
z <- outer(x, y, mf)
z[is.na(z)] <- 4
op <- par(bg = "white")
persp(x, y, z, theta = 60, phi = 20, expand = 0.5)


t=seq(0,2.5,by=0.01)    # Curva (u(t),v(t)) escolhida por acaso.
u<-function(s){s*cos(s)+2}
v<-function(s){-(s^2)*sin(s)/4+1}

#curve(u,0,2.5); curve(v,0,2.5); 
plot(u(t),v(t),'l',col="red")   # curva (u(t),v(t))

persp(x, y, z, theta = 60, phi = 20, expand = 0.5)->res
lines(trans3d(u(t),v(t),mf(u(t),v(t)), res), col = "red", lwd = 2) # Inclui curvas no gráfico em 3d



# Método de Euler oara resolver a equação e'(t)=grad f(e(t)), e(0)=(0.5,1), para determinar numericamente o ponto de máximo.
hh=0.0001
gradmf<-function(u){          # aproximação de f'(x,y) ou do gradiente de f(x,y).
dfx=(mf(u[1]+hh,u[2])-mf(u[1]-hh,u[2]))/(2*hh)
dfy=(mf(u[1],u[2]+hh)-mf(u[1],u[2]-hh))/(2*hh)
c(dfx,dfy)
}

t0=0  # tempo inicial
tf=20  # t final 
e0=c(2.4,0.1) # condição inicial 
n=200
h=(tf-t0)/n  # Tamanho do passo
tt=seq(t0,tf,by=h)
Y=matrix(0,2,length(tt))
Y[,1]=e0
for ( i in 1:(length(tt)-1)){
Y[,i+1]=Y[,i]+h*gradmf(Y[,i])
}


print("Aproximação par aponto de mínimo"); Y[,length(tt)]
print("Aproximação para o valor mínimo"); mf(Y[1,length(tt)],Y[2,length(tt)])

lines(trans3d(Y[1,],Y[2,],mf(Y[1,],Y[2,]), res), col = "blue", lwd = 2) # Inclui curvas

plot(t,mf(u(t),v(t)),'l',col="red")   # curva f(u(t),v(t))

plot(tt,mf(Y[1,],Y[2,]),col = "blue",'l') # curva f(u(t),v(t)) para caminho de maximização obtido pelo método de Euler.


#-----------------------Segunda forma: método de Newton para maximizar e resolver a equação grad f(u)=0.

hh=0.0001
Hessmf<-function(u){          # aproximação de f''(x,y) ou da hessiana de f(x,y).
A=matrix(0,2,2)
v=u; v[1]=u[1]+hh;  w=u;w[1]=u[1]-hh
A[1,]=(gradmf(v)-gradmf(w))/(2*hh)
v=u; v[2]=u[2]+hh;  w=u;w[2]=u[2]-hh
A[2,]=(gradmf(v)-gradmf(w))/(2*hh)
return(A)
}

t0=0  # tempo inicial
tf=20  # t final 
e0=c(0.4,0.1) # condição inicial 
n=200
h=(tf-t0)/n  # Tamanho do passo
tt=seq(t0,tf,by=h)
Y=matrix(0,2,length(tt))
Y[,1]=e0
for ( i in 1:(length(tt)-1)){
w=solve(Hessmf(Y[,i]),-h*gradmf(Y[,i]))
Y[,i+1]=Y[,i]+w
}


print("Aproximação par aponto de mínimo"); Y[,length(tt)]
print("Aproximação para o valor mínimo"); mf(Y[1,length(tt)],Y[2,length(tt)])

plot(tt,mf(Y[1,],Y[2,]),col = "blue",'l') # curva f(u(t),v(t)) para caminho de maximização obtido pelo método de Euler-Newton.