AjusteNaoLinear, R

Language: Layout:

# 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.7543449, 66.1229695, 56.3927738, 47.1209530, 38.9938584, 31.3461193, 
24.2470210, 17.7540897, 11.1669463,  4.0083690, -0.1860567,  3.9043642, 
11.3967455, 17.5330911, 24.1004559, 31.2966078, 38.8839710, 47.1480662, 
56.5175926, 66.3541594, 76.6679233)
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, sendp 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.

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

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

gradg<-function(u){   # função gradiente de g(u).
a=u[1];b=u[2];c=u[3]
p=c(0,0,0)
for (i in 1:n){
p=p+( y[i]- 20*x[i]^2/log(a*x[i]^4+b*x[i]^2+c))*(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)}
p
}

# Teste da função gradg(u)
gradg(c(1,2,3))

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

ZeroEuler<-function(u0,t,m){
w=u0
h=t/m
for (i in 1:m){w=w-h*gradg(w)}
w
}

# Teste de aproximação de solução do problema de minimização. 
u0=c(1,1,2); t=10; m=10000
g(u0)                       # Teste de g(u0).

u=ZeroEuler(u0,t,m) ; u     # Aproximadamente u(t)
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.

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

[ + ] Show input
Absolute running time: 1.22 sec, cpu time: 1.66 sec, memory peak: 47 Mb, absolute service time: 1,43 sec