31-08-2020-Exemplo AjusteNaoL, 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( 0.1617647, 0.1701807, 0.1818182, 0.1985294, 0.2236842, 0.2638889, 0.3333333,
 0.4659091, 0.7500000,  1.3750000, 2.0000000, 1.3750000, 0.7500000, 0.4659091, 0.3333333,
 0.2638889, 0.2236842, 0.1985294, 0.1818182, 0.1701807, 0.1617647)
plot(x,y,col="red")

# Por alguma experiência em lidar com dados, supomos ainda que
# f(s)=(a+bs^2)/(c+ds^2), 
# em que a, b e c são 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(x_i))^2, sendo u=(a,b,c,d).
# 
# 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]; d=u[4]
p=0
for (i in 1:n){
p=p+( y[i]-(a+b*x[i]^2)/(c+d*x[i]^2) )^2}
p/2
}

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

gradng<-function(u){   # função gradiente de g(u).
h=0.00001
p=0*u
m=length(u)
for ( i in 1:m){v=u; v[i]=u[i]+h;p[i]=g(v)-g(u)}
p/h
}

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

#---- 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*gradng(w)}
w
}

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

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

u[1]/u[3]; u[2]/u[3]; u[4]/u[3]  # Comparação comoutro método.

# Definição da aproximação de f(u).

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

curve(fap,-5,5)
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.28 sec, cpu time: 1.75 sec, memory peak: 46 Mb, absolute service time: 1,54 sec