# Gráfico em 3 dimensões

# Resolver a equação f(x,y)=(f1(x,y),f2(x,y))=(0,0).

import math

def f1(u):  # Definição de f1
    p=math.exp(u[1]+u[0])-2
    return p

print('Teste de definição de função f1')
print(f1([1,2]))

def f2(u):   # Definição de f2
    x=u[0]
    y=u[1]
    f2=math.sin(x)+math.cos(y)-1
    return f2
print('Teste de definição de função f2')
print(f2([1,2]))


def f(u):   # Definição de f=(f1,f2)
  p=[f1(u),f2(u)]
  return p
  
print('Teste de definição de função f')
print(f([1,1]))

def crie_matriz(n_linhas, n_colunas, valor):	
	    matriz = [] # lista vazia
	    for i in range(n_linhas):
	        # cria a linha i
	        linha = [] # lista vazia
	        for j in range(n_colunas):
	            linha.append(valor)
	
	        # coloque linha na matriz
	        matriz.append(linha)
	
	    return matriz


def Jacf(u): # Jacobiano de f(x,y)
   x=u[0]
   y=u[1]
   dfx= [math.exp(x+y),math.exp(x+y)]
   dfy= [math.cos(x),-math.sin(y)]
   A=crie_matriz(2,2,0)
   A[0][0]=dfx[0]
   A[0][1]=dfx[1]
   A[1][0]=dfy[0]
   A[1][1]=dfy[1]
   return(A)

print('Teste de definição de função Jacf')
print(Jacf([1,2]))
print(len(Jacf([1,1])))



def iJacf(u): # inversa do Jacobiano  de f(x,y) (caso 2x2)
 A=Jacf(u)
 B=crie_matriz(2,2,0)
 d=A[0][0]*A[1][1]-A[0][1]*A[1][0] # determinante
 B[0][0]=A[1][1]/d
 B[1][0]=-A[1][0]/d
 B[0][1]=-A[0][1]/d
 B[1][1]=A[0][0]/d
 return(B)  

print('Teste de definição de função inversa de Jacf')
print(iJacf([1,1]))

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



e0=[1,1] # para testes
print('condicao inicial')
print('e0',e0)

def prod(A,b):
    C=[]
    for i in range(len(b)):
        p=0
        for j in range(len(b)):
            p += A[i][j]*b[j]
        C.append(p)
    return C   
    
print('Teste de definição  prod')
print(prod(iJacf(e0),f(e0)))

def pnv(t,v):
    p=[]
    for i in range(len(v)):
        p.append(t*v[i])
    return(p)
    
print('Teste de definição  pnv',pnv(2,e0))

def sv(t,v):
    p=[]
    for i in range(len(v)):
        p.append(t[i]+v[i])
    return(p)
    
print('Teste de definição  sv(e0,e0)',sv(e0,e0))
print('e0',e0)

def g(t,u):
    return(sv(f(u),pnv(t-1,f(e0))))
        
print('Teste de definição  g')
print(g(0,e0))

def rungek2(f,u,a,b,n):
    h=(b-a)/n
    r=[]
    p=u
    r.append(u)
    q=f(u)
    print(prod(iJacf(p),q))
    for i in range(0,n):
        print('u_',i,'=',p)
        k1=prod(iJacf(p),q)
        pb=sv(p,pnv(-h,k1))
        k2=prod(iJacf(pb),q)
        p=sv(p,pnv(-h/2,sv(k1,k2)))
        r.append(p)
    return(r)
  
res=rungek2(f,[1,1],0,1,10) 
print('resultado u_',10,'=',res[10][0],res[10][1])

print('teste resultado: f(u_',10,')=',f(res[10]))