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

f<-function(u){A%*%u-z} # f(u)=Au-v

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

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

# Agora é conhecido que o jacobiano de f(u) é A e o gradiente de g(u) é dado por A^*(Au-z).
#
# 
Jacf<-function(u){ A }   # função jacobiano de f(u).


gradg<-function(u){t(Jacf(u))%*%f(u)}
# Teste da função gradg(u)
#gradg(c(1,2,3))

#---- Início da rotin para o método de Euler.


ZeroEuler<-function(u0,m){
Y=u0
for ( i in 1:m){
B=Jacf(Y)
w=t(B)%*%f(Y)
bb=abs( sum( (( t(B)%*%B )%*%w )*w ) )
h=sum(w*w)/bb
if (bb!=0){Y=Y-h*w }
}
Y}

# Teste de aproximação de solução do problema de minimização. 
u0=c(0,0,0);print("Condição inicial");u0; 
m=1000
print("Teste de g(u0)");g(u0)                       # Teste de g(u0).

u=ZeroEuler(u0,m) ; print("Aproximação"); u     # Aproximadamente u(t)
print("Teste de minímo aproximado");g(u)                        # Teste de minímo aproximado

# Definição da aproximação de f(u) linearizada
flinap<-function(s){u[1]*s^4+u[2]*s^2+u[3]} 

curve(flinap,-5,5,ylim=c(0,700))
points(x,z,col="red")        # Teste de ajuste.

plot(x,z-flinap(x),'l')
# Definição da aproximação de f(u).

fap<-function(s){20*s^2/(log(u[1]*s^4+u[2]*s^2+u[3]))} 

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

plot(x,y-fap(x),'l')

# Se quiser pode aumentar t e n ou escolher outra condição inicial 
# para se convencer do resultado obtido.