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

# Por alguma experiência em lidar com dados, supomos ainda que
# f_u(s)=20*s^2/(log(as^4+b2*s^2+c)), 
# em que u=(a, b, c) é o vetor de parâmetros desonhecidos.
# Para determinar os valores dos parâmetros resolvemos o problema de minimização 
#
#        min g(u)
# 
# em que g(u)= sum (y_i-f_u(x_i))^2=sum f_i(u)^2=|f(u)|^2, sendo f_i(u)=y_i-f_u(x_i) e u=(a,b,c).
# 
# Para fazer isso usamos o método de Euler para resolver a equaçõ diferencial 
#
#      u'(t)=-grad g(u(t))
#      u(0)=u_0

# Seguem os comando para as definições de funções necessárias.


f<-function(u){    # função g(u)
	a=u[1];b=u[2];c=u[3]
	p=0*1:n
	for (i in 1:n){	p[i]=y[i]- 20*x[i]^2/log(a*x[i]^4+b*x[i]^2+c) }
	p
}

# Teste do código para f(u).
#f(c(1,2,3)) 

g<-function(u){    # função g(u)
t(f(u))%*%f(u)/2
}

# Teste do código para g(u).
#g(c(1,2,3)) 

Jacf<-function(u){   # função gradiente de g(u).
	a=u[1];b=u[2];c=u[3]
	J=matrix(0,n,3)
	for (i in 1:n){
		J[i,]=(20*x[i]^2/(log(a*x[i]^4+b*x[i]^2+c))^2
		)*( 1/(a*x[i]^4+b*x[i]^2+c) )*c(x[i]^4,x[i]^2,1)}
	J
}
# Teste para o jacobiano de f
#Jacf(c(1,2,3))

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(1,0.5,2.4);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).

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.