# --------------------------------------------------------------------------------------------------
# 3. Resolver o Exercício 2.2.14 das notas de aula.
# Resolver a equação grad gf(x,y,z)=(f1(x,y,z),f2(x,y,z),f3(x,y,z))=(0,0,0), pelo método de Newton  com u_0=(0,0.5,sqrt(2)), em que 
# gf(x,y,z)=y^6-2y^4+y^2+1+2x-6xy^2+3x^2y^4+x^2+2x^3+3x^4y^2+2x^4+x^6+z^4-4z^2+4.

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

f3<-function(u){    # Definição de f3??
    x=u[1]
    y=u[2]
    z=u[3]
    return(4*z^3 - 8*z)
}

#--- Minimo de |grad f(u)|^2 para resolver equação f(u)=0.

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

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

Jf<-function(u){ # A matriz do jacobiano de f(u) é a matriz hessiana de gf(u).
    x=u[1]
    y=u[2]
    z=u[3]
    p=matrix(0,3,3)
    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,0)
    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,0 )
    df3=c(0,0,12*z^2-8)
    p[1,]=df1
    p[2,]=df2
    p[3,]=df3
    p
}

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


# Método de Euler-Newton para a Equação Jf(u(t))u'(t)=-alpha f(u(t))
alpha=9.8

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

# Teste 

print("Resolução equação ")
u0=c(0,0.5,sqrt(2)); t=1; n = 10
u=ZeroEuler(u0,t,n) ; u [,n+1]     # Aproximadamente u(1)
f(u[,n+1])
g(u[,n+1])

# Para classificar o ponto crítico encontrado precisamos olhar para os autovalores 
# da matriz A que é a hessiana da função do exercício, que nesse caso é a matriz do jacobiano que 
# que consta do código.
A=Jf(u[,n+1]);A
# Isso nos diz que um autovalor é 16 (positivo). Os outros autovalores são autovalores da matriz
A[1:2,1:2]
# Como seu determinante é negativo temos que os autovalores também têm sinal oposto.
# Então o ponto crítico encontrado é de sela.
#
# Outra forma de analisar o problema é determinar o polinômio caracteristico da matriz 
# hessiana, por meio das identidades de Newton, como no Exemplo 2.1.12 das notas de aula.
# Isso é feito como 
#Definimos a função que calcula a soma dos elementos da diagonal principal, ou o traço, de uma matriz
tr<-function(A){
p=0; n=length(A[1,])
for ( i in 1:n ){p=p+A[i,i]}
p}	

# O cálculo dos coeficientes do polinômio é feito através  das linhas de comando
	
coef<-function(A){ 
n=length(A[1,]);	a=0*1:(n+1)
s=0*1:n; s[1]=tr(A)
a[1]=1;a[2]=-s[1]; B=A
for ( i in 2:n){B=B%*%A; s[i]=tr(B); q=0
for (j in 1:i){q=q+a[j]*s[i+1-j]}; a[i+1]=-q/i}
a=(-1)^n*a;a}
a=coef(A);a	

# Nesse caso  temos policarac(x)=a[1]x^3 + a[2]x^2+ a[3]x+a[4].
policarac<-function(s){a[1]*s^3 + a[2]*s^2+a[3]*s+a[4]}
curve(policarac,-8,17)	
# Como o polinômio tem raiz negativa, o ponto encontrado é de sela.