#
# Link para o Colab: https://[Log in to view URL]
#
# Resolver a equação f(x,y)=(f1(x,y),f2(x,y))=(0,0) ou o sistema de equações
# f1(x,y)=0
# f2(x,y)=0.
#
# Primeiro vamos definir as funções na linguagem que estamos usando.

f1<-function(u){   # Definição de f1
x=u[1]
y=u[2]
return(2-6*y^2+6*x*y^4+2*x+6*x^2+12*x^3*y^2+8*x^3+6*x^5)
}

f2<-function(u){   # Definição de f2
x=u[1]
y=u[2]
p=6*y^5-8*y^3+2*y-12*x*y+12*x^2*y^3+6*x^4*y
return(p)
}

# Agora vamos fazer os graficos das curvas implícitas f1(x,y)=0 e f2(x,y)=0.

xx <- seq(-2,2, by=0.05)
yy <- seq(-2,2,by=0.05) 
z=matrix(0,length(xx),length(yy))
for ( i in 1:length(xx)){
for ( j in 1:length(yy)){z[i,j]=f1(c(xx[i],yy[j]))}}

contour(xx,yy,z,level=0,col="red")
par(new=TRUE)

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

contour(xx,yy,zz,level=0,col="blue",main="Curvas implicitas")
legend(1.2,1.8, legend=c("f1(u)=0", "f2(u)=0"),
       col=c("red","blue"), lty=c(1,1), cex=0.8)
       

f<-function(u){      # Definição de função f(u), para resolver a equação f(u)=0, 
c( f1(u), f2(u))     # com u=(x,y).
}

# Para a resolução da equação f(u)=f(x,y)=0, vamos usar o valor mínnimo da função 
# g(u)=|f(u)|^2/2.

g<-function(u){
t(f(u)) %*% f(u)/2
} 

# Para a resolução numérica da equação, vamos usar a equação diferencial ordinária
# Jf(u(t))u'(t)=-f(u(t)), com a escolha u(0)=u_0=(1,1), nesse exemplo.
# 
 
Jf<-function(u){  # Jf(u) é matriz do jacobiano de f(u).
x=u[1];y=u[2]
p=matrix(0,2,2)
df1=c(6*y^4+2+12*x+36*x^2*y^2+24*x^2+30*x^4, -12*y+24*x*y^3+24*x^3*y)
df2=c( -12*y+24*x*y^3+24*x^3*y, 30*y^4-24*y^2+2-12*x+36*x^2*y^2+6*x^4 )
p[1,]=df1
p[2,]=df2
p
}

# Método de Euler resolver numéricamente a Equação Jf(u(t))u'(t)=- f(u(t))

ZeroEuler<-function(u0,t,n){
u=matrix(0,2,n+1)
u[,1]=u0
tj=t/n
for (i in 1:n){
h=solve( Jf(u[,i]),-tj*f(u[,i]))    # Resolve o sistema linear Jf(u)h=-tj f(u)
u[,i+1]=u[,i]+h}
u
}

# Teste 

print("Resultado pelo metodo de Newton")
u0=c(1,1); t=6; n = 6            # tamanho de passo tj=1.
u=ZeroEuler(u0,t,n) 
print("Valor da raiz u="); print(u [,n+1])    # Aproximadamente u(6)
print("Valor de f(u)="); print(f(u[,n+1]))
print("Valor de g(u)="); print(g(u[,n+1])[1,1])

xn=u[1,]; yn=u[2,]
plot(xn,yn,'l',col="yellow", main="Euler para Jf(u(t))u'(t)=-f(u(t)), com n=6 e passo h=1")



# Poderíamos resolver a equação f(u)=0 usando o método dos gradientes para g(u)
# Para comparar com o método de Newton, vamos definir o vetor gradiente de g(u).

gradg<-function(u){
p=t(Jf(u))%*%f(u)
t(p)
}

# Método de Euler para a Equação u'(t)=-[Jf(u(t))]*f(u(t)), com u(0)=u_0=(1,1).

ZeroEuler<-function(u0,t,n){
u=matrix(0,2,n+1)
u[,1]=u0
h=t/n
for (i in 1:n){
u[,i+1]=u[,i]-h*gradg(u[,i])}
u
}

# Teste 

print("Resulado pelo metodo dos gradientes ")
u0=c(1,1); t=6; n = 100000        # Note que aqui é necessário que o passo tj seja bem pequeno.
u=ZeroEuler(u0,t,n) 
print("Valor da raiz u="); print(u [,n+1])    # Aproximadamente u(6)
print("Valor de f(u)="); print(f(u[,n+1]))
print("Valor de g(u)="); print(g(u[,n+1])[1,1])


x=u[1,]; y=u[2,]
plot(x,y,col="green",'l',main="Euler para u'(t)=-[Jf(u(t))]*f(u(t)), com n=10000 e passo h=6/n")

# Juntando gráficos

contour(xx,yy,z,level=0,col="red")
par(new=TRUE)

contour(xx,yy,zz,level=0,col="blue",main="Curvas implicitas e trajetorias")
legend(1.2,1.8, legend=c("f1(u)=0", "f2(u)=0"),
       col=c("red","blue"), lty=c(1,1), cex=0.8)
par(new=TRUE)
points(xn,yn,'l',col="yellow", main="Euler para Jf(u(t))u'(t)=-f(u(t)), com n=6 e passo h=1")
points(x,y,'l',col="green", main="Euler para Jf(u(t))u'(t)=-f(u(t)), com n=6 e passo h=1")