# Gráfico em 3 dimensões

# Resolver a equação f(x,y)=(f1(x,y),f2(x,y))=(0,0).

f1<-function(u){   # Definição de f1
r=u[1]
s=u[2]
return((cos(r) + s^3) / sqrt(r))
}

x <- seq(0.1,4, by=0.05)
y <- seq(-4,4,by=0.05)
mf1<-function(r,s){u=c(r,s); f1(u)}# auxiliar para gráfico
z=matrix(0,length(x),length(y))
for ( i in 1:length(x)){
for ( j in 1:length(y)){z[i,j]=f1(c(x[i],y[j]))}}


require(grDevices) # for trans3d
persp(x, y, z, theta = 30, phi = 40, expand = 1)
contour(x,y,z,level=0,col="red")
par(new=TRUE)

f2<-function(u){   # Definição de f2
r=u[1]
s=u[2]
p=tan(s^2 + 3)*cos(r^3)
return(p)
}

zz=matrix(0,length(x),length(y))
for ( i in 1:length(x)){
for ( j in 1:length(y)){zz[i,j]=f2(c(x[i],y[j]))}}

contour(x,y,zz,level=0,col="blue")

legend(2,3.5, legend=c("f1(u)", "f2(u)"),
       col=c("red","blue"), lty=c(1,1), cex=0.8)

persp(x, y, zz, theta = 30, phi = 40, expand = 1)

f<-function(u){c(f1(u),f2(u))} # f(x,y)

mf<-function(u){   # Definição de mf=|f|^2; com coordenadas para gráfico
f1(u)^2+f2(u)^2 
}

print("Teste de mf=|f|^2:"); mf(c(1,2))


gradmf<-function(u){ # gradiente de mf(x,y) aproximado
h=10^(-5)
v=u; v[1]=u[1]+h; v1=u; v1[1]=u[1]-h; 
dfx= mf(v)-mf(v1)
v=u; v[2]=u[2]+h; v1=u; v1[2]=u[2]-h; 
dfy= mf(v)-mf(v1)
p=c(dfx,dfy)/(2*h)
return(p)
} 


print("Teste de gradmf:"); gradmf(c(1,2))


Jacf<-function(u){ # Jacobiano aproximado de f(x,y)
h=10^(-5)

v=u; v[1]=u[1]+h; v1=u; v1[1]=u[1]-h; 
dfx= f(v)-f(v1)

v=u; v[2]=u[2]+h; v1=u; v1[2]=u[2]-h; 
dfy= f(v)-f(v1)

A=matrix(0,2,2)
A[1,]=dfx; A[2,]=dfy
p=A/(2*h)
return(p)
} 

print("Teste do Jacobiano de f:"); Jacf(c(1,2))


#--------------------- Método de Euler para resolver a equação u'(t)=-grad mf(u(t)), e(0)=(2,-2)

t0=0  # tempo inicial
tf=10  # t final 
e0=c(1,-1) # condição inicial 
n=10000
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("Resultados com o Método de Euler para resolver a equação u'(t)=-grad mf(u(t)), e(0)=(2,-2)")
print("|f(u0)|^2"); mf(e0) # teste da escolha
print("Aproximação para o ponto de mínimo ou raiz  de f(u)=0"); Y[,length(tt)]
print("Aproximação para o valor mínimo"); mf(Y[,length(tt)])
print("Aproximação para gradiente de |f|^2no valor mínimo"); gradmf(Y[,length(tt)])
print("Aproximação para Jacobiano de f no valor mínimo"); Jacf(Y[,length(tt)])

#-------------------------------------------------------------

plot(Y[1,],Y[2,],col="blue",'l',main = "Curva solução aproximada para u'(t)=-grad mf(u(t)), e(0)=(2,-2)")

mfY=0*tt
for ( i in 1:length(tt)){mfY[i]=mf(Y[,i])}

plot(tt,mfY,'l', col = "blue",main = "Decaimento de |f|^2")

#------------------------------
# --------------------------- Método de Euler para resolver a equação u'(t)=-[Jacf(u(t))]* x f(u(t)), e(0)=(2,-2)


t0=0  # tempo inicial
tf=10  # t final 
e0=c(1,-1) # condição inicial 
n=10000
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*t(Jacf(Y[,i]))%*%f(Y[,i])
}

print("Resultados obtidos pelo Método de Euler para resolver a equação u'(t)=-[Jacf(u(t))]* x f(u(t)), e(0)=(2,-2)")
print("Aproximação para o ponto de mínimo"); Y[,length(tt)]
print("Aproximação para o valor mínimo"); mf(Y[,length(tt)])

plot(Y[1,],Y[2,],'l',col="green",main = "Curva solução aproximada para u'(t)=-[Jacf(u(t))]* x f(u(t)), e(0)=(2,-2)")

mfY=0*tt
for ( i in 1:length(tt)){mfY[i]=mf(Y[,i])}

plot(tt,mfY, 'l',col = "green",main = "Decaimento de |f|^2")