# Gráfico em 3 dimensões

# Curva
# Superfície z=f(x,y)

x <- seq(-1.2,0.8, by=0.05)
y <- seq(-1.2,1.2,by=0.05)
mf<-function(s,r){z=s+r*1i; Mod(z^4+z^3+z^2+z+2)^2# f(x,y); auxiliar para gráfico
}

require(grDevices) # for trans3d
z <- outer(x, y, mf)
z[is.na(z)] <- 4
op <- par(bg = "white")
persp(x, y, z, theta = 30, phi = 40, expand = 1)
contour(x,y,z,levels=c(0.5,1,2,3,4,5,10))

#
persp(x, y, z, theta = 30, phi = 40, expand = 1)->res

# Método de Euler oara resolver a equação e'(t)=grad f(e(t)), e(0)=(0.5,1)

gradmf<-function(u){ # gradiente aproximado de f(x,y)
s=u[1]
r=u[2]
h=10^(-5)
dnfx= mf(s+h,r)-mf(s-h,r)
dnfy=mf(s,r+h)-mf(s,r-h)
c(dnfx,dnfy)/(2*h)
} 

Hmf<-function(u){ # Hessiana aproximada de f(x,y)
s=u[1]
r=u[2]
h=10^(-5)
v=u; v[1]=u[1]+h; w=u; w[1]=u[1]-h
Hfx= gradmf(v)-gradmf(w)

v=u; v[2]=u[2]+h; w=u; w[2]=u[2]-h
Hfy= gradmf(v)-gradmf(w)
H=matrix(0,2,2)
H[1,]=Hfx; H[2,]=Hfy
return(H/(2*h))
} 


t0=0  # tempo inicial
#tf=10  # t final 
e0=c(4,7) # condição inicial 
mf(e0[1],e0[2]) # teste da escolha
n=16
#h=(tf-t0)/n  # Tamanho do passo


tt=0:n; hh=tt

Y=matrix(0,2,length(tt))
Y[,1]=e0

for ( i in 1:(length(tt)-1)){
w=gradmf(Y[,i])
bb=sum((Hmf(Y[,i])%*%w)*w)
h=sum(w*w)/bb
if (bb!=0){Y[,i+1]=Y[,i]-h*w}
tt[i+1]=tt[i]+h
hh[i]=h
}

print("Ponto de minimo de f(z) ou raiz de p(z)"); Y[,length(tt)]
print("Valor minimo aproximada"); 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(Y[1,],Y[2,],'l',col="blue")


plot(tt,mf(Y[1,],Y[2,]),'l', col = "blue")


plot(tt,hh,'l', col = "green")