# 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
x=u[1]
y=u[2]
f1=y*x^3+y^4+2*x*y+x+y+1 
f1 
}

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

f2<-function(u){   # Definição de f1
x=u[1]
y=u[2]
f2=y*x^2+2*x*y+x+5
f2
}

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

f<-function(u){   # Definição de f=(f1,f2)
c(f1(u),f2(u))
}

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

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))


iJacf<-function(u){ # inversa do Jacobiano  aproximado de f(x,y) (caso 2x2)
A=Jacf(u); B=0*A
d=A[1,1]*A[2,2]-A[1,2]*A[2,1] # determinante
B[1,1]=A[2,2]; B[1,2]=-A[2,1]; B[2,1]=-A[1,2]; B[2,2]=A[1,1]
p=t(B)/d
return(p)
} 

print("Teste da inversa do Jacobiano de f:"); iJacf(c(1,2))

iJacf(c(1,2))%*%Jacf(c(1,2))



alpha<-function(u){min(eigen(t(Jacf(u))%*%Jacf(u))$values)} # valores singulares ao quadrado

print("Teste autovalor do quadrado do valor singular de Jacf"); alpha(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=1  # t final 
e0=c(2,-2) # condição inicial 
n=5000
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)])
print("Aproximação para o inverso do Jacobiano de f no valor mínimo"); iJacf(Y[,length(tt)])

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

plot(Y[1,],Y[2,],col="blue",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, col = "blue",main = "Decaimento de |f|^2",sub = "Estimativa de decaimento de mf como mf(u0) x exp(-beta x t) em vermelho")

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

alpha0=0*tt
for ( i in 1:length(tt)){   # Para teste de decaimento
alpha0[i]=sqrt(alpha(Y[,i]))}
beta=min(alpha0)
print("beta=menor valor singular"); beta

points(tt,mf(Y[,1])*exp(-beta*tt), col = "red",xlim=c(0,1))



plot(tt,sqrt(alpha0), col = "green",main = "Estimativa dos valores singulares de Jacf")


# --------------------------- 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=1  # t final 
e0=c(2,-2) # condição inicial 
n=5000
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,],col="blue",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, col = "blue",main = "Decaimento de |f|^2")



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


t0=0  # tempo inicial
tf=  # t final ????
e0=c(2,-2) # condição inicial 
n=50
      # Tamanho do passo ????


tt=0*1:n; tt[1]=t0

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

for ( i in 1:(length(tt)-1)){
A=Jacf(Y[,i]); B=t(A)%*%A
w=t(A)%*%f(Y[,i])                # w= grad g(u)                     
h= t(w)%*%w /(t(B%*%w) %*%w)     # Escolha de passo ótimo
tt[i+1]=tt[i]+h[1]
Y[,i+1]=Y[,i]-h[1]*t(A)%*%f(Y[,i])
}

print("Resultados obtidos pelo Método de Euler acelerado 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,],col="blue",main = "Curva solução aproximada (acelerada) 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, col = "blue",main = "Decaimento de |f|^2")



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


t0=0  # tempo inicial
tf=2  # t final ????
e0=c(2,-2) # condição inicial 
n=50
h=tf/n   # Tamanho do passo ????
alpha=1
beta=5


tt=seq(t0,tf,by=h)

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

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

for ( i in 1:(length(tt)-1)){
A=Jacf(U[,i]); B=t(A)%*%A
w=t(A)%*%f(U[,i])                # w= grad g(u)                     
V[,i+1]=V[,i]+h*t(A)%*%f(U[,i])-beta*h*V[,i]
U[,i+1]=U[,i]-alpha*h*V[,i+1]
}

print("Resultados obtidos pelo Método de Euler acelerado-b 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"); U[,length(tt)]
print("Aproximação para o valor mínimo"); mf(U[,length(tt)])

plot(U[1,],U[2,],col="blue",main = "Curva solução aproximada (acelerada-b) 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(U[,i])}

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


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

#---------------------------Euler---->Newton



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



print("Resultados obtidos pelo Método de Euler para resolver a equação u'(t)=-[Jacf(u(t))]^(-1) 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,],col="blue",main = "Curva solução aproximada para u'(t)=-[Jacf(u(t))]^(-1) 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, col = "blue",main = "Decaimento de |f|^2")