# 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 
}

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

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

Jacf<-function(u){ # Jacobiano de f(x,y)
x=u[1]
y=u[2]
dfx= exp(x+y)*c(1,1)
dfy= c(cos(x),-sin(y))

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

iJacf<-function(u){ # inversa do Jacobiano  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)
} 

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

#---------------------------Continuacao-homotopia



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



print("Valores a´rpximados por homotopia, com u(0)=(1,1)")
print("u(1) aproximado"); Y[,length(tt)]
print("f(u(1)) aproximado"); f(Y[,length(tt)])

plot(Y[1,],Y[2,],col="blue",main = "u'(t)=-[Jacf(u(t))]^(-1) x f(u(0)), u(0)=(1,1)")

G<-function(s,u){f(u)-f(e0)+s*f(e0)} # Para testar o método


Gg=0*tt
for (i in 1:length(tt)){Gg[i]=sum(abs(G(tt[i],Y[,i])))}
plot(tt,Gg,xlab="t",ylab="G(t,u(t))")

dJacf<-function(u){x=u[1];y=u[2];-exp(x+y)*(cos(x)+sin(y))}
dJ=0*tt
for (i in 1:length(tt)){dJ[i]=dJacf(Y[,i])}
plot(tt,dJ,xlab="t",ylab="det(Jacf(u(t))")