AjusteLinearizado-11-01-21, 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.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")

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

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

# 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.

# Após linearização, tomando z=exp(20*x^2/y), temos que fl_u(s)=as^4+b2*s^2+c 
# é uma família de funções para aproximar os dados linearizados.
#
# Para determinar os valores dos parâmetros resolvemos o problema de minimização 
#
#        min gl(u)
# 
# em que gl(u)= |Au-z|^2, sendo u=(a,b,c) e A a matriz cujas 21 linhas são os vetores (x[i]^4,x[i]^2,1).
# 
# Para fazer isso usamos o método de Euler para resolver a equaçõ diferencial 
#
#      u'(t)=-grad gl(u(t))
#      u(0)=u_0

# Seguem os comando para as definições de funções necessárias.
#
A=matrix(0,21,3)
for (i in 1:n){
A[i,]=c(x[i]^4,x[i]^2,1)}

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

gl<-function(u){t(f(u))%*%f(u)/2}

B=t(A)%*%A; b=t(A)%*%z           # Adaptando o sistema para A*Au=A*v, ou Bu=b, em que B=A*A é positiva definida e b=A*v.

gradgl<-function(u){B%*%u-b}

h=0.0000001                          # Tamanho do passo h, para calcular numéricamante u(20), em que u(0)=(1,1,1,1,1) e u'(t)=-Bu+b.
n=0.25/h
u=c(0.1,0.1,0.1)

for ( j in 1:n){
u=u-h*gradgl(u)                   # Método de Euler
}
print(" Aproximação para u");u                                # Aproximação para u(0.22).
print(" Aproximação para f(u)") ;f(u)                         # Verificação do erro de aproximação.

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

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

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

plot(x,z-fapl(x),'l')

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

# 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.

[ + ] Show input
Absolute running time: 6.13 sec, cpu time: 8.54 sec, memory peak: 57 Mb, absolute service time: 6,36 sec