# Gabarito resumido da Avaliação I

#1. (Exercício 1.4.4 das notas de aula). O exercício pede para usar a aproximação polinomial
# obtida no Exercício 1.4.2 das mesmas notas, para obter uma aproximação para a solução de equação 
# tg(t)-150=0.
# 
# Resolução
# No no Exercício 1.4.2 é dada uma aproximação poliLag(u) para tg(u), apenas para u entre o e 1. Por isso 
# é sugerido o uso da identidade triginométrica tg(1+u)=(tg(1)+tg(u))/(1-tg(1)tg(u)), porque 
# um pouco de noçao de trigonometria nos diz que a solução da equação dada está próxima de pi/2.
pi/2

#-------- Segue os cógigos copiados do link disponível no referido exemplo:
F<-function(s,u){ 1+u^2 } # Condição inicial y(0)=0, y'=1+y^2=F(t,y). 

n=10
a=0; b=1
h=(b-a)/n 

t=seq(a,b,by=h)  # vetor t=(t0, t1,...,tn)
w=0*t 

w[1]=0 # condição incial y(0)=0

for (i in 1:n){ # Runge-Kutta4 para obter aproximações para 10 valores de tg(t_i)
k1=F(t[i],w[i])
k2=F(t[i]+h/3,w[i]+h*k1/3)
k3=F(t[i]+2*h/3,w[i]-h*k1/3+h*k2)
k4=F(t[i]+h,w[i]+h*(k1-k2+k3))
 w[i+1]=w[i]+h*(k1+3*(k2+k3)+k4)/8
}

plot(t,w,col="red")

#------------ Interpolação polinomia
n=length(w)
L<-function(i,s){
p=1
for (j in 1:n){
if ( i!=j){p=p*(s-t[j])/(t[i]-t[j])}
}
p
}

poliLag<-function(s){ # Polinômio na forma de Lagrange
p=0
for ( i in 1:n){p=p+w[i]*L(i,s)}
p
}

curve(poliLag,a,b)
points(t,w,col="red")

#-------------------------------------------------------------------------
# Queremos resolver a equação tg(t)-150=0, mas não podemos usar a 
# função tangente nem sua inversa. Usaremos então a aproximação que obetmos 
# por meio de interpolação e do método de Runge-Kutta.
# Note que o problema proposto pode ser escrito como
#     t(t)=tg(1+u)= (1+tg(u))/(1-tg(u))=150. 
# Usando a aproximação polinomial sugerida escrevemos a equação, de forma aproximada como 
# (poliLag(1)+poliLag(u))/(1-poliLag(1)poliLag(u))-150=0
# 
# Que fica melhor escrita sem fração como 
#
# (poliLag(1)+poliLag(u))-150(1-poliLag(1)poliLag(u))=0
#
# Depois disso a dica é para utilizar conceitos do capítulo das notas de aula que segue o exercício.
# O referido capítuo é sobre máximos e mínimos de funções.
# Por isso tranformamos a equação que acabamos de obter em uma função, como segue
aproeq<-function(u){
(poliLag(1)+poliLag(u))-150*(1-poliLag(1)*poliLag(u))
}

curve(aproeq,0,1,col="red") # raiz aproximada x=0.5


# Isso nos diz que a solução para o problema prosto pode ser vista gráficamente.
# Como o capítulo em questão trata de localização de pontos de mímino e depois usa isso para resolver equações, na Seção 2.2.
#
# Optamos então por resolver a equação dada por meio do método de Newton, usando 
a função que obtemos a pouco.

daproeq<-function(s){
h=10^(-5)
return((aproeq(s+h)-aproeq(s-h))/(2*h))} # Aproximação para a derivada d
# Método de Newton (aproximado)

x=0.5
for ( i in 1:10){x=x-aproeq(x)/daproeq(x)}
x
aproeq(x)
print("Raiz aproximada de  tg(t)-150=0");1+x

#------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------
#
# 2. Exercício 2.1.17 das notas de aula.
# parte a) Um emprestimo de 1200 reais é pago em seis parcelas fixas mensais de 300 reais, sem entrada.
# Qual é a taxa de juros cobrada?
#
# Solução:
# A dívida é um função da taxa de juros cobrada; e deve se anular após seis pagamentos 
# mensais de 300 reais.
# Definimos a função dívida 
divida<-function(r){           # Função dívida 
p=1200                         # valor emprestado
for ( i in 1:6){p=p*(1+r)-300} # juros mensal menos pagamentomensal
return(p) }                    # dívida no último pagamento em função de

curve(divida,0,0.2,col="blue") 
# vemos no gráfico que a dívida será paga após seis parcelas apenas se 
# a taxa de juros estiver próxima de 0.13, ou de 13%.

# Como feito na questão anterior, usamos o método de Newton para determinar uma raiz de 
# da função divida.

r=0.13 # aproximação inicial

dndivida<-function(s){     # derivada numerica
h=10^(-5)
return((divida(s+h)-divida(s-h))/(2*h))} # Aproximação para a derivada de f

# Método de Newton (aproximado)

for ( i in 1:10){r=r-divida(r)/dndivida(r)}

r;divida(r)
# Observe que com essa taxa a dívida fica
d0=1200;do           # Data de empréstimo
d1=d0*(1+r)-300;d1   # 1 mês depois
d2=d1*(1+r)-300;d2   # 2 meses depois
d3=d2*(1+r)-300;d3   # 3 meses depois
d4=d3*(1+r)-300;d4   # 4 meses depois
d5=d4*(1+r)-300;d5   # 5 meses depois
d6=d5*(1+r)-300;d6   # 6 meses depois

# parte b).... Exercício para entregarem até sexta feira.

# --------------------------------------------------------------------------------------------------
# 3. Resolver o Exercício 2.2.14 das notas de aula.
# Resolver a equação f(x,y)=(f1(x,y),f2(x,y),f3(x,y))=(0,0.5,sqrt(2)).

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 |f(u)|^2 para resolver equação f(x)=0.

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

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

Jf<-function(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.