gradconj<-function(A,b,m){ # Método dos gradientes conjugados para o sistema para Au=v, 
                           # com A positiva definida.     
B=A; bb=b           
u=0*b
w=B%*%u-bb                        # w0
v=w                              # d1
alpha= t(w)%*%w /(t(B%*%w) %*%w); alpha=alpha[1]     # Escolha de h0
u=u-alpha*w                       # u1

for ( j in 1:m){
w=(B%*%u-bb)
beta=t(B%*%w) %*%v/(t(B%*%v) %*%v); beta=beta[1]
v=w-beta*v                   # Método de Euler
alpha= t(v) %*%w/(t(B%*%v) %*%v);alpha=alpha[1]
u=u-alpha*v                      # Método quase Euler (integral)
}
return(u)
}

gradconjnp<-function(A,b,m){      # Método dos gradientes conjugados para o sistema para Au=v. 
B=t(A)%*%A; bb=t(A)%*%b           # Adaptando o sistema para A*Au=A*v, ou Bu=b, em que B=A*A é positiva definida e b=A*v.
u=0*b
w=B%*%u-bb                        # w0
v=w                              # d1
alpha= t(w)%*%w /(t(B%*%w) %*%w); alpha=alpha[1]     # Escolha de h0
u=u-alpha*w                       # u1

for ( j in 1:m){
w=(B%*%u-bb)
beta=t(B%*%w) %*%v/(t(B%*%v) %*%v); beta=beta[1]
v=w-beta*v                   # Método de Euler
alpha= t(v) %*%w/(t(B%*%v) %*%v);alpha=alpha[1]
u=u-alpha*v                      # Método quase Euler (integral)
}
return(u)
}

# Dados do exercício
h=0.25
x0=-0.75
x=seq(x0,1,by=h);print("x");x
y=c(9.649763, 10.807633, 12.502805, 15.000000, 18.702063, 24.225860, 32.521764, 45.062646);print("y");y

plot(x,y)
#--------------------
# Resolução seguindo dica que acompanhava o exercício
#
# Resolução do sistema linear para determinar lambda
n=length(x);n

A=matrix(0,4,4)

for (i in 1:4){
A[i,]=y[i:(i+3)]
}
print("primeiro sistema linear");A

v=y[5:8];v

lambda=gradconj(A,-v,8);print("lambda");lambda
#--Teste de precisão de lambda como solução do sistema linear
A%*%lambda+v

#--Definição do polinômio com coeficientes lambda
polino<-function(s)
{
lambda[1]+lambda[2]*s+lambda[3]*s^2+lambda[4]*s^3+s^4
}

# Derivada do polinômio
dpolino<-function(s)
{
lambda[2]+2*lambda[3]*s+3*lambda[4]*s^2+4*s^3
}

curve(polino,0.75,1.75)

curve(dpolino,0.75,1.75)

# Método de Newton
beta=c(0.7,1.1,1.4,1.7) # Condição inicial para o método de Newton, obtida graficamente.
for ( i in 1:100){beta=beta-polino(beta)/dpolino(beta)}
print("beta");beta

#--Teste de beta como raiz do polinômio
polino(beta)

# Matriz para determinar alpha
B=matrix(0,4,4)
for (i in 1:4){
for (j in 1:4){
B[i,j]=beta[j]^(i-1)}}
print("segundo sistema linear");B;y[1:4]

alpha=gradconjnp(B,y[1:4],10);print("alpha");alpha

#--Tesde de precisão de alpha como solução do sistema linear
B%*%alpha-y[1:4]

# Cáclulo de a e de b.
b=log(beta)/h;print("b");b
a=alpha*(beta)^3;print("a");a

# Definição da função interpoladora encontrada
fp<-function(s){a[1]*exp(b[1]*s)+a[2]*exp(b[2]*s)+a[3]*exp(b[3]*s)+a[4]*exp(b[4]*s)}

#--Teste de ajuste da função encontrada com os dados
curve(fp,-1,1)
points(x,y,col="red")

#--Teste de escala do erro cometido, nos pontos de interpolação
plot(x,y-fp(x),col="red")