# Ajuste de curvas por mínimos quadrados.

# Dados os vetores x e y, supomos que y_i=f(x_i)+erro_i, como erro_i um erro "pequeno".
x=seq(-5,5,by=0.5);x  
n=length(x)
y=c(76.8, 66.3, 56.4, 47.2, 38.9, 31.4, 24.4, 17.84, 11.1, 3.7, -0.2, 4.1, 11.21, 17.6,	
24.1, 31.3, 39.0, 47.33, 56.53, 66.13, 76.73)
plot(x,y,col="red")

# Mudança de variável z=e^(20x^2/y)=ax^4+bx^2+c

z=exp(20*x^2/y);z

plot(x,z,col="red")

# Agora cálculos mostram que a função g(u)= |Au-v|^2 / 2
# em que A é a matriz e v=z é um vetor, como montamos no que segue.

A=matrix(0,n,3); v=0*x;
for ( i in 1:n){A[i,]=c(x[i]^4,x[i]^2,1);v[i]=z[i]}

gmv<-function(u){  #  |Au-v|^2 / 2
w=A%*%u-v; t(w)%*%w/2}

gmv(c(1,6,3))     # Teste de g(u), com mudança de variável.

# Agora é conhecido que o gradiente de g(u) é dado por A^*(Au-v).
#
# Como o problema é linear, vamos resolver pelo Método dos gradientes conjugados 
B=t(A)%*%A; B; b=t(A)%*%v;b      # Adaptando o sistema para A*Au=A*v, ou Bu=b, 
# em que B=A*A é positiva definida e b=A*v.

n=10
u=c(2,2,2)                       # u0
w=B%*%u-b                        # w0=grag(u0)
d=w                              # d1
h= t(w)%*%w /(t(B%*%w) %*%w)     # Escolha de h0
u=u-h[1]*w                       # u1
w=B%*%u-b  


for ( j in 1:(n-1)){
a= t(B%*%d) %*%w/(t(B%*%d) %*%d) # Escolha de alpha
d=w-a[1]*d                       # w - projeção ortogonal de w na direção de d
bt= t(d) %*%w/(t(B%*%d) %*%d)    # Escolha de beta
u=u-bt[1]*d                       # (quase) Método de Euler u_{n+1}=u_n-h (w_n-proj_d(w_n))
w=B%*%u-b                        # gradg(u_n)   
}
u                                # Aproximação para u(t), sendo t a soma dos passos h.
A%*%u-v 

gmv(u)

fapmv<-function(s){u[1]*s^4+u[2]*s^2+u[3]}  # Aproximação da função linearizada.

curve(fapmv,-5,5)                           # Teste de ajuste
points(x,z,col="red")


fap<-function(s){20*s^2/(log(u[1]*s^4+u[2]*s^2+u[3]))}  # Voltando a variável original

curve(fap,-5,5,ylim=c(0,78))                            # Teste de ajuste
points(x,y,col="red")