# 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=exp(x+y)-2
f1 
}

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

f2<-function(u){   # Definição de f1
x=u[1]
y=u[2]
f2=sin(x)+cos(y)-1
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)/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 acelerado-b 


t0=0  # tempo inicial
tf=10  # t final ????
e0=c(1,1) # condição inicial 
n=10000
h=tf/n   # Tamanho do passo ????
alpha=10
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)){                  
V[,i+1]=V[,i]+h*t(Jacf(U[,i]))%*%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)=-alphav, v'(t)=[Jacf(u(t))]* x f(u(t))-beta x v(t), e(0)=(0,0), v(0)=(0,0)")
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",'l',main = "Curva solução aproximada (acelerada-b) para u(t)")

plot(V[1,],V[2,],col="blue",'l',main = "Curva solução aproximada (acelerada-b) para v(t)")

mfY=0*tt;Mfb=0*mfY
for ( i in 1:length(tt)){mfY[i]=mf(U[,i])
Mfb[i]=mf(U[,i])+alpha*(V[1,i]^2+V[2,i]^2)/2
}

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


plot(tt,Mfb, col = "blue",'l',main = "Decaimento de |f|^2+alpha x |v|/2")